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FOREWORD 

The  present  report  consists  of  two  volumes  which  describe  an  approxi- 
mate nonlinear  analysis  of  solid  rocket  motors  and  T-burners  and  the 
associated  computer  programs.  Volume  I contains  the  analytical  basis  for 
the  computer  programs  and  the  results  of  the  parametric  studies,  while 
Volume  II  describes  the  computer  programs  and  serves  as  a user's  manual. 

The  Investigation  Is  entitled  APPROXIMATE  NONLINEAR  ANALYSIS  OF  SOLID 
ROCKET  MOTORS  AND  T-BURNERS.  The  two  volumes  are  additionally  subtitled 
as  follows: 

Volume  I - Analysis  and  Results 
Volume  II  - Computer  Program  User's  Manual 

This  Investigation  was  sponsored  by  the  Air  Force  Rocket  Propulsion 
Laboratory,  Edwards  AFB,  California  93523  under  contract  number  F04611-75-C-0036 
with  Capt.  Jack  Donn  as  technical  monitor.  Program  management  was  provided 
by  B.  T.  Zlnn  (Co-principal  Investigator)  while  project  engineering  was 
provided  by  E.  A.  Powell  (Co-principal  Investigator)  and  M.  S.  Padmanabhan 
(Post  Doctoral  Fellow). 

This  report  has  been  reviewed  by  the  Information  office  (01)  and  is 
releasable  to  the  National  Technical  Information  Service  (NTIS) . At  NTIS, 
it  will  be  available  to  the  general  public.  Including  foreign  nations. 

This  technical  report  has  been  reviewed  and  is  approved  for  publication. 


WILLIAM  F.  MORRIS,  Colonel,  USAF 
Chief,  Technology  Division 
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1. 


INTRODUCTION 


Experience  with  solid  propellant  rocket  motors  indicates  that  the 
full-thrust  operation  of  such  motors  is  never  completely  steady.  A certain 
degree  of  rough  combustion  (i.e.,  low  amplitude  random  flow  oscillations)  is 

usually  present  which  does  not  seriously  affect  system  performance.  On 

occasion,  larger  amplitude  organized  oscillations  may  develop  in  the  com- 
bustion chamber  and  nozzle  flow  field  during  the  ignition  transient  or  after 

full  thrust  has  been  obtained^.  Such  organized  oscillations  are  evidence  of 

a positive  feedback  between  the  unsteady  propellant  combustion  process  and 
the  flow  field  disturbance.  These  instabilities  may  be  classified  as  either 
spontaneous  instabilities,  which  originate  from  flow  or  combustion  noise,  or 
pulsed  instabilities,  which  result  from  the  Introduction  of  a sufficiently 
large  disturbance  in  a motor  which  is  stable  with  respect  to  small  amplitude 
disturbances.  It  is  observed  that,  after  passing  through  a transient  period, 
both  spontaneous  and  pulsed  instabilities  reach  a limiting  amplitude  (or 
limit-cycle)  at  which  they  oscillate  with  a frequency  that  is  close  to  the 
frequency  of  one  of  the  chamber's  acoustic  modes.  These  limit-cycle  oscilla- 
tions are  often  characterized  by  nonsinusoidal  pressure  waveforms  that  have 
sharp  peaks  and  flattened  minima.  The  existence  of  limit-cycles , pulsed 
instability,  and  nonsinusoidal  waveforms  are  all  caused  by  the  nonlinearities 
of  the  system. 

Small  amplitude  combustion  instability  may  occur  without  detrimental 
effects  to  the  system  while  large  amplitude  oscillations  may  lead  to  increase 
in  mean  chamber  pressure  and  burning  rate,  excessive  heat  transfer  rates,  and 
severe  vibration  levels.  The  occurrence  of  any  one  of  these  may  result  in 
malfunction  or  destruction  of  the  rocket  motor.  Thus,  it  is  important  for 
the  rocket  designer  to  be  able  to  estimate  the  limiting  amplitude  of  the 
pressure  oscillations  as  well  as  the  conditions  under  which  pulsed  instability 
may  occur.  To  do  this,  a theoretical  approach  capable  of  determining  the 
nonlinear  stability  characteristics  of  solid  propellant  rocket  motors  is 
required. 


^Dehority,  G.  L.  and  Price,  E.  W.,  "Axial  Mode,  Intermediate  Frequency 
Combustion  Instability  in  Solid  Propellant  Rockets,"  Naval  Weapons  Center 
Tech.  Publ.  5654,  October  1974. 
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Reliable  theoretical  approaches  capable  of  predicting  the  character- 
istics of  combustion  instabilities  and  the  conditions  under  which  they  are 
most  likely  to  occur  are  badly  needed.  To  be  of  practical  use,  these 
analyses  should  be  conceptually  simple  and  easily  adaptable  for  use  by 
engine  designers.  In  addition,  these  techniques  should  be  capable  of  solv- 
ing combustion  instability  problems  without  exceeding  memory  core  limita- 
tions of  current  computers  and  without  requiring  excessive  ciimputation  time. 
These  considerations  are  particularly  relevant  to  multi -dimensional  com- 
bustion instability  problems.  As  the  above  requirements  cannot  be  satisfied 
by  any  of  the  "exact"  numerical  solution  techniques,  (e.g.,  finite-difference 
or  finite-element  methods)  one  must  resort  to  the  use  of  one  or  more  appro- 
ximate solution  techniques.  The  development  of  such  solution  techniques  for 
the  analysis  of  nonlinear  axial  combustion  instabilities  in  solid  rockets  is 
the  subject  of  this  report. 

In  recent  years,  several  nonlinear  combustion  instability  theories 

2 3 4 5 

have  been  developed.  Kooker  and  Zinn  ’ and  Levine  and  Culick  ’ used 
finite-difference  techniques  to  obtain  numerical  solutions  to  the  conserva- 
tion equations  describing  the  two-phase  flow  of  gases  and  particles  in  solid 
rocket  combustors.  However,  limitations  on  available  computer  core  size  and 
computational  time  make  it  impractical  to  use  these  so-called  "exact" 
analyses  in  parametric  studies  and  engine  design.  Hence,  efforts  have  been 
undertaken  in  recent  years  to  develop  reliable  approximate  solution  techni- 
ques that  are  free  from  the  above-mentioned  limitations.  In  the  first 
known  attempt  at  the  development  of  such  an  approximate  technique,  Zinn, 
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Kooker,  D.  E.  and  Zinn,  B.  T.,  "Numerical  Solution  of  Axial  Insta- 
bilities in  Solid  Propellant  Rocket  Motors,"  Proceedings  of  the  10th  JANNAF 
Combustion  Meeting.  CPIA  Publication  243,  Vol.  I,  December  1973,  pp . 389-415. 
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Kooker,  D.  E.  and  Zinn,  B.  T. , "Numerical  Investigation  of  Nonlinear 
Axial  Instabilities  in  Solid  Rocket  Motors,"  BRL  CR  141,  March  1974. 
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Levine,  J.  N.  and  Culick,  F.E.C.,  "Numerical  Analysis  of  Nonlinear 
Longitudinal  Combustion  Instability  in  Metalized  Propellant  Solid  Rocket 
Motors,"  Proceedings  of  the  9th  JANNAF  Combustion  Meeting.  CPIA  Publication 
231,  Vol.  1,  December  1972,  pp.  141-163. 

^Levine,  J.  N.  and  Culick,  F.E.C.,  "Nonlinear  Analysis  of  Solid  Rocket 
Combustion  Instability,  AFRPL-TR-74-45,  Vol.  I,  October  1974. 
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Powell  and  Lores  have  used  the  Galerkin  method  in  the  analysis  of  nonlinear 

combustion  Instabilities  in  liquid  rockets  with  quasi-steady  nozzles  exper- 

6 7 8 

iencing  longitudinal  and  transverse  instabilities.  Powell  and  Zinn 

later  extended  these  theories  to  situations  in  which  the  instabilities  are 

three-dimensional  and  the  rocket  combustors  are  attached  to  conventional 

nozzles.  More  recently,  a solution  procedure  utilizing  the  Galerkin  method 

9 10 

together  with  the  method  of  averaging  has  been  applied  by  Culick  ’ in  the 
analysis  of  nonlinear  instabilities  in  solid  rockets.  While  the  above-men- 
tioned approximate  theories  appear  promising,  more  work  is  needed  to  deter- 
mine their  reliability  and  range  of  applicability,  since  these  theories  are 
limited  by  the  assumptions  upon  which  they  are  based. 

The  main  objective  of  this  investigation  is  to  evaluate  the  usefulness 

of  approximate  nonlinear  analysis  techniques  for  predicting  the  stability 

behavior  of  solid  rocket  motors.  This  objective  has  been  accomplished 

through  a comparison  of  the  predictions  of  the  approximate  model,  based  on 

the  methodology  developed  by  Zinn,  Powell,  Lores^  ^ and  Culick^’ with  the 

2 3 4 5 

more  exact  analysis  developed  by  Kooker  and  Zinn  ’ and  Levine  and  Culick  ’ 
Therefore,  e,i.Lensive  development  of  nonlinear  models,  such  as  developing  new 
theory  or  programming  a large  analysis,  was  not  necessary.  However,  modifi- 
cations to  both  the  approximate  and  "exact"  analyses  were  made  in  order  to 
accomplish  the  objectives  of  this  program.  These  modifications  will  be 
discussed  in  Sections  2 and  3 of  this  report. 


^Loras,  M.  E.  and  Zinn,  B.  T.,  "The  Prediction  of  Nonlinear  Longitudinal 
Combustion  Instability  in  Liquid  Propellant  Rockets,"  NASA  CR-120904, 

April  1972. 

^Zinn,  B.  T.  and  Powell,  E.  A.,  "Nonlinear  Combustion  Instability  in 
Liquid  Propellant  Rocket  Engines,"  Proceedings  of  the  13th  Symposium 
(International)  on  Combustion,  The  Combustion  Institute,  1971,  pp.  491-503. 
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Powell,  E.  A.  and  Zinn,  B.  T. , "The  Prediction  of  Nonlinear  Three- 
Dimensional  Combustion  Instability  in  Liquid  Rockets  with  Conventional 
Nozzles,"  NASA  CR-121279,  October  1973. 
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Culick,  F.E.C.,  "Nonlinear  Behavior  of  Acoustic  Waves  in  Combustion 
Chambers,"  Proceedings  of  the  10th  JANNAF  Combustion  Meeting.  CPIA  Publica- 
tion 243,  Vol.  I,  December  1973,  pp.  417-436. 

Culick,  F.E.C.,  "Nonlinear  Behavior  of  Acoustic  Waves  in  Combustion 
Chambers,"  California  Institute  of  Technology  Report,  April  1975. 
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The  Investigation  is  divided  into  four  tasks.  In  the  first  task  the 
influence  of  the  gasdynamic  nonlinearities,  which  result  in  mode-coupling, 
was  assessed.  For  this  task  all  processes  except  gasdynamical  mode-coupling, 
were  described  by  linear  models.  In  the  second  task  the  approximate  and 
"exact"  analyses  were  extended  to  include  the  effects  of  nonlinear  particle 
damping  and  nonlinear  combustion  driving.  In  the  third  task  the  approximate 
analysis  was  applied  to  determine  the  nonlinear  stability  characteristics  of 
T-burners,  Finally,  the  restrictions  and  limitations  on  the  application  of 
the  method  of  averaging  to  the  solution  of  solid  rocket  instability  problems 
was  investigated  in  the  fourth  task. 

The  development  of  the  approximate  analysis  is  given  in  Section  2, 
while  the  modifications  to  the  "exact"  analysis  are  discussed  in  Section  3. 

The  results  obtained  during  the  performance  of  the  four  tasks  described 
above  are  presented  and  discussed  in  Section  4.  These  results  include  an 
extensive  parametric  study  to  assess  the  Importance  of  gasdynamic  mode 
coupling,  particle  size  and  concentration,  propellant  response  function,  and 
particle  and  combustion  nonlinearities  upon  growth  and  decay  rates,  frequencies, 
limiting  amplitudes  and  waveforms  for  motors  and  T-burners.  In  many  cases, 
the  approximate  solutions  are  compared  with  the  corresponding  "exact"  solu- 
tions and  with  experimental  data.  Conclusions  drawn  from  the  results  of 
this  investigation  are  presented  in  Section  5,  along  with  recommendations 
for  improvement  of  the  approximate  model. 
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2.  APPROXIMATE  ANALYSIS 

2. 1 Basic  Assumptions  and  Conservation  Equations 

The  solid  rocket  configuration  investigated  herein  is  described  in 
Figure  1.  A full-length  cylindrical ly-perforated  solid  propellant  grain 
containing  aluminum  particles  burns  unsteadily  inside  a cylindrical  combus- 
tor-nozzle combination.  The  combustion  process  takes  place  inside  a thin 
region  immediately  adjacent  to  the  propellant.  The  products  of  combustion 
are  assumed  to  consist  of  a mixture  of  a single  gaseous  species  and  aluminum 
oxide  particles  which  enter  the  engine  core  flow  at  the  outer  edge  of  the 
combustion  zone  and  acquire  axial  velocity  through  Interaction  with  the 
engine  core  flow.  During  instability,  the  interaction  between  the  wave 
motion  and  the  combustion  process  results  in  unsteady  burning  which  pro- 
vides the  energy  needed  to  sustain  the  instability.  At  the  same  time  the 
interaction  between  the  wave  motion  and  the  particles,  the  core  flow,  and 
the  nozzle  results  in  wave  energy  dissipation.  Limit-cycle  conditions  are 
achieved  when  the  wave  energy  addition  over  a cycle  is  balanced  by  the  wave 
energy  removal  over  a cycle. 

The  following  assumptions  are  used  in  the  development  of  the  theoreti- 
cal model:  (1)  the  flow  in  the  engine  is  one-dimensional  and  it  consists  of 
a single  gaseous  species  and  spherical  particles  that  can  be  characterized 
by  a single  average  disimeter;  (2)  the  gas  phase  is  thermally  and  calorically 
perfect;  (3)  the  gas  phase  is  inviscid  and  non-heat  conducting;  (4)  the 
thermal  energy  transfer  between  gas  and  particle  phases  is  neglected;  (5) 
the  momentum  exchange  between  the  gas  phase  and  the  particles  due  to  viscous 
interaction  can  be  described  by  Stokes'  Drag  Law^^;  (6)  the  burning  rate 
responds  to  pressure  oscillations  only,  and  it  can  be  described  by  a linear 
burning  response  function;  (7)  the  interaction  of  the  oscillations  in  the 
combustor  with  the  flow  in  the  nozzle  can  be  adequately  described  by  using 

g 

an  appropriate  nozzle  admittance  ; (8)  the  Mach  number  of  the  mean  flow  is 
small;  and  (9)  only  moderate -amplitude  disturbances  are  considered  in  this 
analysis . 

Assumptions  (1),  (2),  (3),  (7),  (8)  and  (9)  are  commonly  used  in 


^^Landau,  L.  D.  and  Lifshitz,  E.  M.,  Fluid  Mechanics,  p.  66,  Pergamon 
Press,  1959. 
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combustion  instability  analyses^  Assumption  (4)  could  be  relaxed,  as 

was  done  in  Reference  10,  at  the  expense  of  other  assumptions,  but  it  is 
introduced  here  to  simplify  the  analysis.  The  use  of  Stokes'  Law  and  linear 
combustion  response  (assumptions  (5)  and  (6))  are  consistent  with  the  require- 
ment that  only  the  effect  of  gasdynamical  nonlinearities  is  considered  in  the 
analysis  of  nonlinear  mode-coupling  (first  task).  The  effect  of  particle 
and  combustion  nonlinearities  upon  the  resulting  instability  are  considered  in 
the  second  task,  and  the  results  are  presented  later  in  the  report. 

Using  assumptions  (1)  through  (4)  above,  the  system  of  nondimens ional 

conservation  equations  that  describe  the  unsteady  behavior  of  the  combustor 

3 

two-phase  flow  can  be  expressed  in  the  following  form  : 

Mass  Conservation  (gas  and  particulate  phases) 


5p  3 2 m 

5^  + ^ (pu)  = 


(1) 


3 p 3 2 m 

T?  + 3l^  (Pp  = ~1^ 


(2) 


Momentum  Conservation  (gas  and  particulate  phases) 
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Equation  of  State  (gas  phase) 


P = P. 


(6) 


In  the  above  equations  the  quantity  describes  the  momentum 

exchange  between  the  gas  phase  and  the  particulate  matter.  For  linear 
particle  damping,  it  is  derived  from  the  Stokes  Drag  Law  for  a single 
sphere  in  a steady  flow;  this  F^  is  given  by 


F = KP  (u  - u) 
P P P 


(7) 


where 


18m. 

* *2 
pa  c 
'^m  o 


(8) 


The  case  of  nonlinear  particle  damping  is  discussed  later  in  this  section. 
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The  unsteady  mass  burning  rate  is  described  by  the  following  linear 
12 

response  function  law  (pressure  coupling  only): 


g 


(9) 


where  is  the  complex  response  function  that  must  be  provided  from  either 

experimental  data  or  an  analysis  of  the  unsteady  combustion  process.  In 

this  report  the  response  function  k?  is  taken  to  be  the  well-known,  two- 

12 

parameter  (i.e.,A,B)  expression  given  by 


R = 


n A B 

X +f  - (1+A)  + AB 

A. 


(10) 
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where  A and  B are  related  to  the  steady-state  combustion  of  the  propellant 
n is  the  pressure  exponent  in  the  steady-state  burning  law  and  \ is  a complex 
solution  of  the  quadratic  equation 


X = ) - i n 


(11) 


where  Q is  a nondimens ional  frequency  parameter.  A typical  curve  of  real 
part  of  the  response  function  versus  frequency  parameter  i)  is  shown  in 
Figure  2.  Nonlinear  expressions  for  unsteady  mass  burning  rate  are  given 
later  in  this  section. 

2 . 2 Order  of  Magnitude  Analysis  and  Approximate  Wave  Equations 

To  proceed  with  the  analysis,  each  dependent  variable  is  expressed  as 
a sum  of  a space-dependent  steady-state  quantity  and  a perturbation  which  is 
both  time  and  space  dependent.  Following  the  methodology  introduced  in 

l2 

Culick,  F.E.C.,  "A  Review  of  Calculations  for  Unsteady  Burning  of  a Solid 
Propellant,  "AIAA  J. . Vol,  6,  No.  12,  pp.  2241-2244,  December,  1968. 
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References  (6),  (7)  and  (8)  each  perturbation  quantity  and  the  mean  flow  Mach 

number  are  assumed  to  be  of  0(e)  where  e is  a small  ordering  parameter  that  is 

a measure  of  the  wave  amplitude.  Substituting  the  assumed  expressions  for  the 

dependent  variables  into  Equations  (1)  through  (6),  neglecting  all  terms  of 
2 

order  higher  than  e and  subtracting  out  the  steady-state  equations  result 
in  the  derivation  of  the  desired  second  order  perturbation  equations.  Intro- 
ducing the  velocity  potentials  $ and  defined  by 


u = r — 
p dx 


into  the  resulting  system  of  second  order  wave  equations,  using  the  substitu- 
13  14 

tion  principle  ’ to  eliminate  variables,  the  system  of  second  order  equations 
can  be  combined  to  yield  the  following  nonlinear  wave  equation  describing  the 
unsteady  two-phase  flow  inside  the  motor: 


^xx  - ^t  = 2 ^ ^xt  + (Y+1)  ^ ^ + (v-1) 

+ - ^%>t]  - °p  "p  ^%>xt  - "^c  ^ \ (^3^ 

Considering  the  nature  of  the  problem  under  consideration,  an  in- 
spection of  the  above  equation  suggests  that  it  has  the  proper  form.  The 
left-hand  side  of  the  equation  is  the  wave  equation  operator;  the  first  two 
terms  on  the  right-hand  side  describe  the  effects  of  the  mean  flow;  the  third 
and  fourth  terms  describe  gasdynamical  nonlinearities  while  the  fifth  and 
sixth  terms  describe  the  gas-particle  interaction;  finally,  the  seventh  term 
represents  the  "driving". provided  by  the  unsteady  combustion  process. 


^Lighthill,  M.  J.,  Surveys  in  Mechanics,  p.  250,  1956. 

^^Blackstock,  D.  T. , "Approximate  Equations  Governing  Finite-Amplitude  Sound 
in  Thermo-Viscous  Fluids",  General  Dynamics  GD/E  Report  GD-1463-52,  1963. 
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A detailed  derivation  of  Equations  (13)  and  (14)  is  given  in  Appendix  A. 

In  the  analysis  for  a solid  rocket  motor,  the  head  end  (x=0)  of  the 
chamber  is  assumed  to  consist  of  a rigid  wall  requiring  zero  velocities  for 
both  the  gas  and  particles  at  this  location.  The  influence  of  the  exhaust 
nozzle  on  the  wave  motion  in  the  combustion  chamber  is  introduced  through  the 
boundary  condition  at  x=l,  which  is  the  location  of  the  nozzle  entrance  plane. 

g 

As  is  customary  in  combustion  instability  analyses  , the  nozzle  effect  is  in- 
troduced through  an  appropriate  nozzle  admittance  relation.  Thus,  the 
solutions  to  Equations  (13)  and  (14)  are  required  to  satisfy  the  following 
boundary  conditions. 


(15a) 


$ + =0  X = 1 

X ' t 


(15b) 


where  Y = Y + iY.  = (u^/p')  is  the  nozzle  admittance  coefficient.  The 
r 1 x=l 

boundary  conditions  at  the  ends  of  a T-burner  are  described  later  in  this 
section. 

By  making  an  additional  approximation,  the  wave  equation  for  the  gas 
phase  (i.e..  Equation  (13)  ) and  the  particle  potential  equation  (i.e..  Equation  ^ 

J 
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L 


(14)  ) can  be  combined  to  obtain  a single  equation  describing  the  gas  oscilla- 
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tions.  This  approach  is  similar  to  that  employed  by  Culick  in  which  the 

spatial  derivatives  of  the  particle  properties  are  neglected  in  the  particle 

momentum  equation.  The  resulting  first  order  linear  equation  for  u'  is 

P 

solved  analytically  to  obtain 


Ke 


-Kt 


u (x,t  )e  dt 


(16) 


where  short-term  transient  effects  have  been  neglected.  This  solution  is 
substituted  into  the  gas  phase  equation,  eliminating  the  need  to  numerically 
solve  an  additional  equation  describing  the  particle  oscillations  (i.e., 
Equation  (14)  ).  The  resulting  wave  equation  for  the  gas  phase  then  becomes: 

i - $ = Zu  i + (Y+l)  4^  + 2$  f + (\-l)?>  $ 

XX  tt  xt  dx  t X xt  t XX 


+ KD  - (V-l)Kp  u $ - K 0 ? 

P t ' p X ‘"p  ^ 


^ ,,3-  -Kt  r ^ Kt'  , 
+ K p e $ e dt 

D J 


^ „2-  - -Kt  r , Kt'  , 

+ (\-l)  K 0 u e I f e dt 

P W X 


- v«h  4^  4 

c dx  t 


(17) 


A detailed  derivation  of  Equation  (17)  is  also  given  in  Appendix  A. 

Most  of  the  results  presented  in  this  report  were  obtained  by  an 
approximate  solution  of  Equations  (13)  and  (14)  for  which  fewer  assumptions 
were  necessary.  Approximate  solutions  of  Equation  (17),  however,  are  also 
presented  in  order  to  determine  the  conditions  under  wlnich  solving  the  single 


13 


equation  yields  accurate  results,  since  a considerable  saving  in  computer 
time  is  expected  by  using  this  approach. 


2.3  Application  of  the  Galerkin  Method 

The  approximate  solutions  of  Equations  (13)  and  (14)  or  Equation  (17) 
are  obtained  by  means  of  the  Galerkin  method,  which  is  a special  case  of  the 
method  of  weighted  residuals^^’ Before  proceeding  with  the  application  of 
the  Galerkin  method  to  this  problem,  a brief  description  of  the  method  will  be 
given. 

According  to  the  "classical"  Galerkin  method  the  dependent  variables 

are  expanded  in  terms  of  a set  of  functions  that  identically  satisfy  the 

imposed  boundary  conditions.  The  proper  choic-  of  the  functions  Y^  is  aided 

by  information  trom  such  sources  as  experimental  data,  the  solution  of  a 

linearized  version  of  the  same  problem,  or  frcxn  solutions  of  closely  related 

problems.  Each  of  the  Y is  multiplied  by  an  arbitrary  constant  or  function 

n 

yet  to  be  determined.  These  expansions  are  then  substituted  into  the  governing 
differential  equations  to  form  residuals,  and  the  arbitrary  constants  (or 
functions)  are  determined  by  imposing  the  condition  that  the  residuals  of  the 
differential  equations  be  orthogonal  to  all  the  functions  Y^^ . 

When  the  boundary  conditions  which  the  solution  must  satisfy  are  compli- 
cated, it  is  usually  impossible  to  find  a set  of  functions,  Y , that  can 

n 17,18 

identically  satisfy  these  boundary  conditions.  It  has  been  shown  that 
when  this  occurs  it  is  possible  to  obtain  approximate  solutions  by  properly 
combining  the  boundary  residuals  with  the  differential  equations'  residuals 
when  applying  the  orthogonality  conditions.  Thus,  the  resulting  orthogonality 
conditions  have  the  following  form: 


^^Finlayson,  B.  A.  and  Scriven,  L.  E.,  "The  Method  of  Weighted  Residuals-- 
A Review,"  Applied  Mechanics  Reviews,  Vol.  19,  No.  9,  September  1966,  pp. 
735-744. 

^^Ames,  W.  F.,  Nonlinear  Partial  Differential  Equations  in  Engineering, 
Academic  Press,  1965. 

^^Zinn,  B.  T.  and  Powell,  E.  A.,  "Application  of  the  Galerkin  Method  in  the 
Solution  of  Combustion  Instability  Problems,"  Proceedings  of  the  19th 
International  Astronautical  Congress,  Vol.  3,  1970,  pp.  59-73.  (Also 
appeared  in  lAF  Paper  No.  P69.  Originally  published  in  1968.) 

1 Q 

Zinn,  B.  T.  and  Powell,  E.  A.,  "The  Galerkin  Method  and  its  Use  in  the 
Solution  of  Combustion  Instability  Problems",  Proceedings  of  the  5th  ICRPG 
Combustion  Conference,  December  1968,  pp.  138-144. 
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(18) 


J 


V 


(T.v.'p)  dV 


0 i = 1 ,2,3 


where  E and  B.  respectively  represent  the  differential  equation  and  boundary 
i 1 

residuals.  The  residuals  are  obtained  when  the  approximate  expansions  lor 
the  density,  velocity  vector,  and  pressure  are  substituted  into  the  equations 
and  boundary  conditions  expressing  conservation  of  mass,  momentum,  and  energy. 
The  above  method  of  solution  accounts  for  the  presence  of  mass,  momentum,  and 
energy  sources  both  within  the  flow  field  and  on  the  boundaries  of  the  rocket 
combustor. 

In  order  to  obtain  approximate  solutions  of  Equations  (13)  and  (14) 
by  the  Galerkin  method,  the  velocity  potentials  are  expressed  as  series  of 
acoustic  modes  with  unknown  time-dependent  amplitudes  as  shown  below: 

N 

? = Xi  A (t)  X (x) 

j=l  J J (19a) 


N 

T = X B (t)  X (X) 

P J J 


(19b) 


where  N is  the  number  of  modes  included  in  the  series  expansion  and  X^  (x) 
are  the  acoustic  eigenfunctions  for  the  axial  modes.  For  a rocket  motor, 
the  X 's  are  the  axial  acoustic  modes  for  the  chamber  with  no  mean  flow,  a 

j 

rigid  wall  boundary  condition  at  the  head-end,  and  a nozzle  admittance 
boundary  condition  at  the  other  end.  These  functions  (complex)  are  given  by: 


X.  (X) 


cosh 


(ibjX) 


(20) 
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where  are  the  corresponding  axial  acoustic  eigenvalues  (ctimplex).  The 

unknown  time-dependent  amplitude  functions  (i.e.,  A.(t)  and  B (t)  ) are 

J j 

complex  functions  of  time.  The  assumed  form  for  the  chosen  expansions 
has  been  guided  by  related  work  in  this  area  and  experimental  evidence^ ’ ^ ^ 
showing  that  rocket  instabilities  usually  involve  the  acoustic  modes  of  the 
combustor. 

The  assumed  expansions  are  substituted  into  Equations  (13)  and  (14) 

and  boundary  condition  (15b)  to  form  the  residuals  E {$,$  ] , E ],  and 

P P P 

B {i);  any  of  these  residuals  becomes  identically  zero  when  the  assumed 
solutions  identically  satisfy  the  corresponding  differential  equations  or 
boundary  conditions.  To  determine  the  unknown  time-dependent  coefficients, 
these  residuals  are  required  to  satisfy  the  following  Galerkin  orthogonality 

J.  ■ 1 

conditions  : 


1 

I E I Xj  (x)  dx  - B I $ I (1)  = 0 

o 


(21a) 


1 

J ^ “j  <’■>  " 

o 


j = 1.  2. 


N 


(21b) 


where  X (x)  are  the  cimiplex  conjugates  of  the  axial  acoustic  eigen- 

j 

functions  in  Equations  (14).  Performing  the  integrations  indicated  in  Equation 
(21)  results  in  a system  of  coupled  second-order  nonlinear  ordinary  differentia 
equations  describing  the  behavior  of  the  time-dependent  amplitudes.  These 
equations  are  similar  in  form  to  those  developed  in  References  6 through  10, 
and  they  can  be  expressed  as  follows: 


2N 


d^A 


, dt 


m) 


ITF=  1 


N N 


-l  dA  V’  r 

+ "5?  + L L l^iCj.'n.n) 


m=  1 n=  1 


+ D2( 

j,n>,n) 

A 

m 

* 

dA 
n 

dt 

+ D^(j,m,n) 

* 

A 

m dt 

+ D^( 

j,m,n) 

•k 

A 

m 

■jt 

dA 
n 

dt 

II 

o 

= 1,2. 

(22) 


, dA  , 

L dT  ^ V = ° j = ^>2>  •••  •^ 

nr=  1 


(23) 


In  the  above  equations,  the  B.(t)  functions  for  j = 1,2,  . . . N have  been 
denoted  by  A^(t)  for  j = N + 1,  . . . 2N;  there  are  2N  unknown  time  dependent 
functions  A.(t)  and  2N  equations  (N  of  Equation  (22)  and  N of  Equation 
(23)). 


The  nonlinear  terms  involving  the  complex  conjugates  appearing 
in  Equation  (22)  arise  because  the  wave  equation  (i.e..  Equation  (13)) 
must  be  modified  for  use  with  the  assumed  complex  solution  given  by  Equations 
(19).  This  modification  is  necessary  because  only  the  real  part  of  the 
assumed  solution  is  physically  meaningful.  It  can  easily  be  shown  that 
ifj  = cp  + i'i'isa  solution  to  Equation  (13),  the  real  part,<p  , is  not  a 
solution  to  Equation  (13).  This  failure  of  cp  to  satisfy  Eq.  (13)  is  due 
to  the  presence  of  the  nonlinear  terms  in  this  equation.  It  can  also  be 
shown,  however,  that  a modified  wave  equation  can  be  constructed  for  which 
the  real  part  of  its  solution  satisfies  the  original  wave  equation  (i.e.. 
Equation  (13).  This  modified  wave  equation  is  derived  in  Appendix  B and 
is  given  by  Equation  (B-8).  Thus  Equation  (22)  was  actually  derived  by 
applying  the  Galerkin  method  to  Equation  (B-8)  rather  than  Equation  (13). 
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The  coefficients  of  the  various  linear  and  nonlinear  terms  appearing 
in  Equations  (22)  and  (23)  arise  from  the  spatial  integrations  indicated 
in  Equations  (21).  Some  of  these  coefficients  are  functions  of  the  steady 
state  properties  u(x),  u^(x),  and  Pp(x)  which  are  obtained  by  solving  the 
steady  state  conservation  equations.  All  of  the  coefficients  involve 
integrals  of  products  of  two  or  three  axial  acoustic  eigenfunctions  with 
various  weighting  functions,  integrated  over  the  length  of  the  chamber. 
Expressions  for  these  coefficients  are  given  in  Appendix  C,  and  the  steady 
state  solutions  are  presented  in  Appendix  D. 

In  deriving  Equation  (23)  by  applying  the  Galerkin  method  to  Equa- 

2 

tion  (14),  the  nonlinear  term  l/2($^)^  was  dropped  in  order  to  be  consistent 
with  the  assumption  of  linear  particle  damping.  As  a result  of  neglecting 
this  term.  Equation  (23)  is  a linear  equation. 

To  determine  the  stability  characteristics  of  a solid  rocket  motor, 
the  form  of  the  initial  disturbance  is  specified  and  the  subsequent  behavior 
of  the  individual  modes  is  determined  by  numerically  integrating  Equation 
(22)  and  (23).  Once  the  time-dependence  of  the  individual  modes  (i.e., 
Aj(t)'s)  is  known,  the  gas  and  particle  potentials,  $ and  4^,  are  cal- 
culated from  Equations  (19).  The  pressure  perturbation  at  any  axial  location 
within  the  chamber  is  related  to  the  real  parts  of  § and  (i,e.,  cp 

and  qip)  by  the  following  approximate  momentum  equation  (Appendix  A); 


CPj.  + u(> 


12  12  d 

X)  cp^  + 2 - 2 'Pt  + KOp(x)  (cp-cPp)  + ^ cp_ 


The  Galerkin  method  was  also  used  to  obtain  approximate  solutions 

of  Equation  (17)  which  was  derived  by  neglecting  the  spatial  derivatives 

of  the  particle  velocity  in  Equation  (14).  The  resulting  differential 

equations  governing  the  behavior  of  the  unknown  amplitude  functions, 

A.(t)  is  given  in  this  case  as  follows: 

J 

- r 1 

1 ^ ( j ,m)  A^  + 02(3  ,m)  + h R C,  (j,m)  + 

ni=l  ' Hr  L ^ ^ J 


+ C^(j,m)  e 


-Kt 


N N 


dA 


dA 


2^  L ( D^a.m.n)  A^ 

m=l  n=l 


+ D,(j.m,n)  A^  ^ + D^(j.m.n)  A^  JT  ) = ° 


dA 
r 

‘m  dt 


(25) 


j = 1,2,  ...  N 


which  is  a system  of  N equations  in  N unknowns.  The  coefficients 


, C^,  and  differ  from  those  in  Equations  (22)  and  are  given  in  Appen- 


dix C;  the  coefficients  D^,  , D^,  and  are  the  same  as  those  appearing  in 

Equations  (22).  In  order  to  handle  the  integral  terms  in  Equations  (25), 


auxiliary  variables  G (t)  are  defined  by 

m 


G (t)  = f A (t')e’^‘^  dt' 

m j in 


(26) 


The  G^(t)'s  are  described  by  the  following  equations 


dG 


— = A^(t)  - KG^(t)  m = 1,2,  ...  N 


(27) 


which  are  obtained  by  differentiating  Equations  (26).  Solutions  for  the 


Aj(t)'s  are  then  obtained  by  numerically  solving  Equations  (25)  simultaneously 


with  Equations  (27).  The  pressure  perturbation  is  then  obtained  from  the 
appropriate  form  of  the  momentum  equation  (Appendix  A): 


P = - Y 


[cp 


+ u(x)cp 


12  1 2 ^ dG 

+ ^ CP„  - n CP  + — 


j.  • ■ 2 ” 2 '^t  dx 


Kp 


^(x)cp  - K^0p(x)e  J cp  e*^*" 


(28) 
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1-  •r-r'C 


where  9 is  the  real  part  of  $ given  by  Equation  (19a)  and  the  integral 
terra  is  given  by; 

t N 

e'^*^  j 'Pe'^'^dt  = Re  G^(t)  X^(x)j  (29) 

o m=  1 


2 . A . Method  of  Averaging 

The  theory  presented  thus  far  represents  a two-stage  siraplif ication 
of  the  original  problem.  In  the  first  stage  the  problem  has  been  reduced 
to  the  solution  of  a pair  of  nonlinear,  partial  differential  equations 
(i.e.,  Equations  (13)  and  (14)).  In  the  second  stage  the  solution  was 
expanded  in  a series  of  acoustic  modes  with  time-dependent  amplitudes  and 
the  Galerkin  Method  was  used  to  replace  the  solution  of  tlie  nonlinear 
partial  differential  equations  with  the  solution  of  a system  of  nonlinear 
ordinary  differential  equations  (i.e.,  Equations  (22)  and  (23)  or  (25)), 

As  demonstrated  in  Reference  10  a further  simplification  of  the  problem  (and 
reduction  in  computation  time)  may  be  obtained  by  using  the  Method  of  Averag- 
ing (referred  to  hereafter  as  MOA)  to  reduce  Equations  (22)  and  (23)  or 
Equations  (25)  into  a system  of  first  order,  coupled  differential  equations 
describing  the  growth  or  decay  of  the  unknown  amplitudes  Aj(t).  A brief 
description  of  the  application  of  the  MOA  to  the  present  problem  is  given 
below;  further  details  are  available  in  References  10  and  19. 

To  apply  the  MOA  the  governing  equations  must  be  expressed  in  the 
following  form; 


2. 

+ ju  . A . 
J J 


F. 

J 


j = 1,  2, 


(30) 


dA. 

where  F.  is  a nonlinear  function  of  the  various  A.  and  -5—^  . 

J J dt 

to  the  MOA  each  unknown  amplitude  A^(t)  is  expressed  as 
Aj(t)  = 3j(t)  sin  (uUjt)  + bj(t)  cos 


According 


(31) 


where  a.(t)  and  b.(t)  are  assumed  to  be  slowly-varying  functions  of  time 
J J 

(i.e.,  their  fractional  changes  during  one  period  are  small).  The  functions 
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cos  (ul;^  t ^ )dt  ^ 


t+Ti 

irV  J sin(x  t')dt' 

j 1 , 


j = 1,2,  ...  N 

where  t ^ = 2n/u)^ 

The  solutions  of  Equations  (32)  require  less  computation  time  as 

these  equations  only  describe  the  slowly-varying  parts  of  the  unknown 

amplitudes  Aj(t),  as  the  oscillatory,  rapidly-varying  parts  of  the  solutions 

are  specified  by  sin  oj^t  and  cos  uj^t  in  Equation  (31). 

The  MOA  was  first  applied  to  Equations  (22)  and  (23),  but  the  resulting 

solutions  were  unstable  when  particles  were  present  in  the  combustor  flow. 

Thus  Culick's  approach^^  was  adopted  and  the  MOA  was  applied  to  Equations 

(25).  These  equations  are  a system  of  coupled  differential  equations  to 

be  solved  for  the  unknown  complex  amplitude  functions,  A^(t).  In  order 

to  apply  the  MOA  to  this  system,  they  must  first  be  separated  into  their 

real  and  imaginary  parts.  This  is  done  by  assuming  that  ~ Ej(t) 

+ iG.(t),  substituting  into  Equations  (25),  and  separating  real  and  imaginary 

parts  to  obtain  the  equivalent  system  of  real  differential  equations  that 

describe  the  behavior  of  the  F.'s  and  G.'s.  Since  these  equations  contain 

J J 

twice  as  many  unknown  functions  (i.e.,  F.  and  G ^ ) as  Equations  (25)  it 
is  convenient  to  re-index  the  unknown  functions  and  their  coefficients 
as  follows: 

F (t)  = B,  ,(t) 

P 2p-l 

G (t)  = B„  (t)  (33) 

P 2p 

Thus  the  B's  with  odd  indices  correspond  to  the  real  parts,  F^(t),  and 
the  B's  with  even  indices  correspond  to  the  imaginary  parts,  G^(t).  The 
corresponding  set  of  differential  equations  is  of  the  following  form: 


1 


> • • • 


) 


(34) 


dB^  dB2 
®2N»  dc~  ’ dT” 


j = 1,2,...  2N 

The  above  equations  are  coupled  in  the  second  derivatives;  that 

I 

is,  there  are  two  or  more  Cq  terms  in  each  equation.  This  coupling  results 
from  the  non-orthogonality  of  the  axial  eigenfunctions.  In  order  to  apply 
the  MOA  to  Equations  (34),  they  must  be  decoupled  by  transforming  to  the 


I £s . / dB^  dB2 

= fj^B2,B2,...  B2^,  ^ 


This  transformation  is  accomplished  by  a matrix  inversion  technique  described 
in  Appendix  C.  Finally  Equations  (35)  can  be  written  in  the  form: 


2 


where 


,m)  B^  + €2(3, m)  + (j,m) 


+ hc^jj^  E2  (j,m)  — +C2(j,m)e 


2N  2N 


D(j,m,n)B^j^ 


j = 1,  2,  ...  2N 


Thus  Equations  (25)  have  been  transformed  to  the  required  form  given  by 
Equations  (30).  The  real  coefficients  C^,  , C3,  E^,  E2.  and  D are  derived 

from  the  complex  coefficients  appearing  in  Equation  (25)  during  the  process 
of  separating  real  and  imaginary  parts  and  decoupling  the  second  derivatives 
(Appendix  C). 

The  MOA  is  then  applied  to  Equations  (36)  and  (37)  by  assuming 

that 


® ^ ~ sin  (oi.t)  + h.(t)  cos  (o.'.t) 

J J J J 

where  the  g^(t)  and  h^(t)  are  slowly-varying  functions  of  time.  In  per- 
forming the  integrations  indicated  in  Equations  (32)  the  g.  and  h . are 
assumed  constant  (the  basic  assumption  is  that  g.  and  h.  do  not  vary  sig 
nificantly  during  the  interval  of  averaging).  The  resulting  equations 
describing  the  time-variation  of  the  functions  g^  and  h^  are  given  by: 


Ifj 

dt 


2N 

2 . 

uu . h . - / h 6 . 

m=i 

Cj  (j  ,m)  + 2 2 '"3*^j 

K + (jj 

m . 

2N 

r_. 

X.  a'  g 6 

m m pi 

m«l 


C2  (j,m)  + h^  (J,m)  + h^  E^  ( j ,m)  - 


K + u) 


-Kt  - M 

•2  (e  ^ - 1) 


' 1 r 1 

— 2 - K cos(ou.t)  — 

K + (T  L J J J J n 


2N  C^Cj.m) 

£ [V"'"’ ■"'’-H 


2N  2N 

S X]  ^ (j.m.n) 
m=l  n=l 


n 

~ ®m  ^n  B2(-t.P.'l)  + — *^m®n  ^1 


(39a) 
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dh^ 

dt 


1 i 2 

■ 2»j  ) “1  *J  ■ ""  ‘pi 


K + LU_ 


2N 

+ uu  h 6 . 

m m p-C, 

m=l 


Cj(j ,m) 


C2(j.m)  + h^R^E^(j,m)  + h^R^E^Cj.m)  - j 


K + 0) 


m 


+ e 


-Kt 


r 2nK  -|  r 1 “^1 

e - 1 K 3in(uL) . t)  + ciu . cos(uu.t)  — 

L 1 JL  J J 


2 2 
K +u). 
J 


2N 


g C3(j,m) 


m=l 


3"-^*  ^ r 

„2  2 [' 

K.  + ou  '- 


- Kh  (0) 


2N  2N  r u)  cu  -j  I 

n=l  n=l  > 


j = 1,2, ...2N 


(39b) 

(39b) 


In  these  equations  6p^  is  zero  for  p and  unity  for  p = <t  and  0^ 

and  0^  are  integrals  of  products  of  three  trigonometric  functions  as  defined 
in  Appendix  C.  The  indices  t,  P>  q are  the  axial  mode  numbers  corresponding 
to  the  indices,  j,  m,  n respectively.  This  notation  is  necessary  because 

two  B.(t)  functions  are  needed  to  describe  each  mode. 

J 

To  obtain  nonlinear  solutions  for  a rocket  motor,  initial  amplitudes 
gj(0)  and  hj(0)  are  specified  and  Equations  (39)  are  integrated  numerically. 
The  resulting  solutions  for  gj(t)  and  hj(t)  yield  the  growth  or  decay  rates 
of  the  various  modes  directly.  To  obtain  the  pressure  waveforms  Equation 
(38)  is  used  to  compute  the  corresponding  which  are  then  combined 

to  obtain  :p  . Once  qj(x,t)  is  known  the  pressure  perturbation,  p'  , is  cal- 
culated using  Equation  (28). 

A few  comments  regarding  computation  time  will  now  be  made.  Apply- 
ing the  Galerkin  method  (without  using  the  MOA)  to  Equation  (17)  yields 
a set  of  2N  second  order  equations  and  2N  first  order  equations  after  the 
real  and  imaginary  parts  have  been  separated.  This  is  equivalent  to  a 
set  of  6N  first  order  equations.  On  the  basis  of  the  number  of  equations 
alone  the  MOA  yields  only  a modest  decrease  in  computation  time.  The  advan- 
tage of  the  MOA  lies  in  the  fact  that  the  g^  and  h^  functions  are  slowly 
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the  numerical  integrations  resulting  in  a substantial  saving  in  computation 
time.  In  order  to  evaluate  the  accuracy  of  the  MOA,  solutions  obtained 
by  the  Galerkin  method  (without  averaging)  and  solutions  obtained  by  the 
MOA  will  be  compared  in  Section  IV. 


2.5  Application  to  Motors 


The  approximate  analyses  (Galerkin  method  and  MOA)  may  be  applied 
to  a variety  of  motor  and  T-burner  configurations.  The  application  of 
the  approximate  methods  to  determine  the  nonlinear  stability  behavior  of 
a motor  with  a full-length  tubular  propellant  grain  and  a quasi-steady 
nozzle  (Figure  1)  is  discussed  in  this  subsection.  The  application  to 
T-burners  is  covered  in  the  next  subsection. 

As  mentioned  previously,  the  coefficients  that  appear  in  Equations 
(22)  and  (23),  (25),  and  (39)  are  dependent  upon  the  steady-state  proper- 
ties u(x),  u (x),  and  p (x).  Therefore,  the  steady-state  solutions  must 
P P 

be  obta^.ied  before  the  approximate  analysis  can  be  applied.  Although  these 

steady-state  quantities  could  be  obtained  by  numerically  integrating  the 

steady-state  versions  of  Equations  (1)  through  (6),  it  is  more  convenient 

to  use  an  approximate  analytical  solution  for  the  steady-state  properties. 

In  deriving  the  approximate  unsteady  equations,  it  was  assumed  that  the 

steady-state  values  of  p , p,  and  h (or  T)  were  constants  (the  variation 

- 2 

in  these  properties  is  of  order  M^  ) . To  this  order  of  approximation  the 
steady-state  gas  and  particle  velocities,  u(x)  and  u^(x),  are  linear  func- 
tions of  the  axial  coordinate  x,  and  the  particle  density  p (x)  is  a con- 
stant. The  following  expressions  were  determined  by  satisfying  the  steady- 
state  continuity  and  momentum  equations  with  p(x)  = 1 and  dp^/dx  = 0: 


u(x)  = 


X = u X 
e 


Up(x)  = 


1 +T/1  + 


(42) 


nip 

■1 

/ 

1 — _■ 

j l+-y  1+  8u  /K 

1+  ’ 

■{  ^ 

C 
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The  derivation  of  these  equations  is  given  in  Appendix  D. 

In  order  to  obtain  approximate  stability  solutions  for  motors, 
the  nozzle  admittance  Y must  also  be  specified.  For  conventional  nozzles 
Y is  a complex  number  which  is  dependent  on  nozzle  geometry  and  oscillation 
frequency.  For  short  nozzles  the  unsteady  nozzle  flow  is  nearly  quasi- 
steady and  Y becomes  a real  number  which  is  independent  of  frequency.  The 
value  of  Y for  a quasi-steady  nozzle  is  determined  by  the  condition  that 
the  Mach  number  at  the  nozzle  entrance  is  a constant.  Assuming  small  per- 
turbations and  neglecting  the  effect  of  particles  yields  the  following 
expression  for  Y: 


Y 

r 


u 

e 


Y.  = 0 
1 


(43) 


2.6  Application  to  T-Burners 

The  approximate  methods  described  in  this  section  were  also  used 
to  analyze  the  nonlinear  behavior  of  T-burners.  Since  the  T-burner  geometry 
differs  in  several  important  respects  from  that  of  motors,  a detailed  dis- 
cussion of  these  differences  is  given  in  this  subsection. 

The  T-burner  geometry  considered  in  this  investigation  is  shown 

* 

in  Figure  3.  The  T-burner  consists  of  a cylindrical  tube  of  length  L with 
★ 

a vent  of  length  at  the  center.  In  the  simplest  configuration,  propellant 
disks  are  placed  at  each  end  of  the  T-burner  (end-burning  onlyl.  For  metal- 
lized propellants  it  is  often  necessary  to  include  short  tubular  grains 
★ 

of  length  L at  the  ends  in  order  to  increase  the  burning  propellant  surface 
° * 
area  (cup  grains).  The  cross-sectional  area  of  the  burner  is  S , while 

the  cross-sectional  area  of  the  cylindrical  grain  is  S As  in  the  case 

c • 
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* 

of  motors,  Lho  axial  variable  x is  normalized  by  the  chamber  length, 

★ * ... 
thus  X = X /L  and  the  dimensionless  lengths  of  the  cylindrical  grains 

and  the  center  vent  are  given  by 


* * * * 
e = / L , / L 


(44) 


For  convenience  in  the  analysis,  the  T-burner  is  divided  into  five 
regions  as  shown  in  Figure  3.  In  Regions  1 and  5 at  the  ends  of  the  burner, 
the  flux  of  burned  gases  and  particles  at  the  lateral  propellant  surface 
is  represented  by  mass  and  energy  source  terms  in  the  conservation  equations. 
Thus  these  regions  are  similar  to  the  interior  of  a motor,  except  for  the 
boundary  fluxes  at  the  ends  (x  = 0 and  x = 1).  In  Regions  2 and  4 there 
is  no  combustion  at  the  lateral  boundary,  so  the  corresponding  source  terms 
are  zero.  Finally  in  Region  3,  the  mass  flux  through  the  center  vent  is 
represented  by  mass  and  energy  sinks  in  the  governing  equations.  For  the 
case  of  end  burning  disks  only,  ^ = 0 and  Regions  1 and  5 vanish. 

In  order  ,o  apply  the  approximate  analysis  techniques  to  T-burners, 
the  analysis  must  be  modified  to  account  for  the  following;  (1)  the  presence 
of  burning  propellant  at  the  ends  of  the  chamber,  (2)  the  presence  of  partial 
length  lateral  propellant  surfaces,  and  (3)  the  presence  of  the  center 
vent.  These  points  will  now  be  discussed  separately. 

End  Burning.  The  burning  propellant  grains  at  the  ends  of  the 
T-burner  must  be  treated  as  boundary  conditions  rather  than  source  terms 
as  in  the  case  of  lateral  grains.  These  burning  surface  boundary  conditions 
replace  the  rigid  wall  boundary  condition  (head-end)  and  nozzle  admittance 
boundary  condition  used  in  the  motor  analysis.  The  boundary  conditions 
at  the  ends  of  the  T-burner  can  be  described  most  conveniently  in  terms 
of  the  response  function  defined  by  Equation  (9): 

-H  YUj^(ft-  l)^j.  =0  at  X = 0 

atx  = l (45) 

where  Uj^  is  the  steady-state  velocity  of  tlie  burned  gases  leaving  the  propellant 
surfaces.  In  deriving  Equation  (45)  from  Equation  (9)  the  assumption  that 
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D = I and  p = 1 was  used,  and  p'  was  replac'?d  by  its  first  order  equi- 
valent The  difference  in  the  signs  for  and  Bj^  arises  because 

a positive  pressure  perturbation  ( < 0)  yields  a positive  velocity  per- 
turbation at  the  left  end  01  '’‘'O  negative  velocity  perturbation 

at  the  right  end  0).  In  other  words,  a positive  flux  (i.e.,  leaving 

the  surface)  gives  positive  velocities  at  x = 0 and  negative  velocities 
at  X = 1,  since  the  positive  direction  for  velocities  is  taken  to  he  the 
direction  of  increasing  x. 

In  applying  the  Galerkin  method  to  the  T-burner,  a boundary  residual 
arises  at  both  ends  of  the  chamber,  rather  than  just  at  one  end  as  in  the 
motor  case.  Thus  the  Galerkin  orthogonality  conditions  given  by  Equations 
(21)  are  replaced  by  the  following  expressions: 


I ~ ~ * \~l  * 

E ’ X.  (x)dx  + B § ' X.  (o) 

/ P\  J J 


B '?^x'^(l)  = 0 

^1  I J 


\ I -k 

E X. (x)dx  = 0 

Jo  P / P 1 J 


(46a) 


(46b) 


j — 

where  the  E and  E are  the  residuals  of  Equations  (13)  and  (14)  and  the 
P 

boundary  residuals  B^  and  B^  are  givo:i  by  Equations  (45).  The  eigenfunc- 
tions X.j  are  still  given  by  Equation  (20),  but  the  corresponding  b^'s  are 
diffi?rent  from  those  used  previously  for  the  analysis  of  motors.  The  pro- 
per selection  of  the  b.'s  of  course,  should  be  the  acoustic  eigenvalues 
for  a chamber  with  boundary  conditions  given  by  Equations  (45).  In  the 
case  of  the  T-burner,  however,  the  velocities  of  the  gases  leaving  the 
propellant  surfaces  ai  ■ sufficiently  small  that  these  eigenvalues  can  be 
approximated  by  the  eigenvalues  for  a closed-ended  chamber.  Thus  the  eigen- 
values b.  are  given  by 


bj  = j n j = l,2,i,...,N  (47) 

Partial  Length  Cylindrical  Grains.  The  lateral  burning  surfaces 
are  treated  as  source  terms  as  was  done  previously  for  motors,  but  the 
approximate  analysis  has  been  extended  to  handle  partial  length  grains. 

Two  grain  configurations  are  considered  in  this  study:  (1)  disk  grains  only 
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•md  (2)  cup  Brains  (Figure  3).  In  the  latter  case  flush  grains  (S  = S ) , 

CO  c 

ire  assumed.  Variable  cross-sectional  area  has  an  important  eftect  on  mode  - 

shape  and  frequency.^'*  I'hus,  to  accommodate  protruding  or  recessed  grains,  1 

the  mode  shapes  used  in  the  assume<l  series  expansion  must  be  those  given 

20  • ■ i 

by  Derr,  which  are  much  more  complicated  than  those  presently  used. 

n 

Thus  a major  modification  of  the  present  analysis  would  be  necessary  to 

accommodate  T-burners  of  variable  cross-sectional  area,  which  is  beyond  ; 

the  scope  of  this  project.  Thus  the  approximate  T-burner  analysis  is  re-  j 

stricted  to  flush  grain  configurations.  _] 

In  applying  the  Galerkin  method  to  analyze  T-burners  with  cup  grains,  3 

the  spatial  integrations  indicated  in  Equations  (46)  must  be  performed 

over  three  distinct  types  of  regions.  In  the  first  type  of  region,  the  * 

source  terms  present  in  the  wave  equation  are  due  to  lateral  surface  burn- 
ing (Regions  1 and  5)  and  Equation  (13)  is  used  to  obtain  the  residual  ^ 

E(^).  In  the  second  type  of  region,  the  source  terms  are  zero  (Regions  ^ 

2 and  4),  and  the  residual  E(§)  is  obtained  by  dropping  the  terms  in 
Equations  (13).  The  third  type  of  region,  the  vent  region,  will  be  discussed 
later.  These  differences  between  the  first  two  types  of  regions  are  auto- 
matically taken  into  account  by  introducing  the  appropriate  steady-state 
solutions  when  the  integrals  are  performed. 

The  steady-state  solutions  for  Regions  1 and  2 (as  well  as  Regions 
4 and  5)  were  obtained  from  the  steady-state  continuity  and  momentum  equa- 
tions. In  order  to  obtain  analytical  expresions  for  u(x),  u^(x),  and 
p (x)  a number  of  simplifying  assumptions  were  made.  These  assumptions 
are  consistent  with  the  order  of  magnitude  approximations  used  in  deriving 
Equations  (13)  and  (14).  A discussion  of  these  approximations  and  the 
derivation  of  the  steady-state  relations  are  given  in  Appendix  D. 

In  the  regions  of  lateral  surface  burning,  the  steady-state  solutions 
are  similar  to  those  for  motors,  except  for  differences  caused  by  the  finite 
velocity  at  the  ends  of  the  burner.  In  Region  1 (0  < x < 0/2)  the  steady- 
state  solutions  are: 


20 

Derr,  R.  L. , "Evaluation  of  a Variable  Area  T-Burner  for  Metallized  Pro- 
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u(x)  = u^(l  ) 


(48) 


Up(x) 


TH-2 


(49) 


Pp(x)  = 


1 + ~n 


1 + T1 


1 + 


2 -,  - mi\ 

rxJ  ti  I 


(50) 


where  Tj  = 4u^/RK.  For  x = 0 these  relations  yield  u(0)  = Uj^ , Up(0)  = , 

and  Op(0)  = C^,  which  agree  with  the  boundary  values  imposed  on  the  steady 
state  solutions.  The  gas  velocity  is  seen  to  vary  linearly  with  x with  a 
slope  given  by: 

dx  R (51) 


while  the  particle  velocity  approaches  a linear  variation  with  x after 
a short  transitional  distance  6x.  Similarly,  the  particle  density  approaches 
a constant  value  tor  x > 6x.  These  asymptotic  values  are  given  as  follows; 


Up(x) 


1+Tl 


(52) 


p^(x)  = (1+H)C^ 


(53) 


For  particle  sizes  of  a few  microns,  T]  is  a small  positive  number  and  the 

exponents  appearing  in  Equations  C49)  and  (50)  are  large,  thus  Equations 

(52)  and  (53)  are  very  good  approximations  for  all  but  very  small  values 

of  X.  in  Region  5,  the  directions  of  the  velocities  u(x)  and  u (x)  are 

P 

opposite  from  those  in  Region  1,  and  the  magnitudes  of  the  steady-state 
properties  are  obtained  from  equations  (48)  through  (50)  by  replacing 
X with  1-x.  For  1 ~8/2<x<l  - 6x  Equations  (52)  and  (53)  are  also 
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valid,  with  u(x)  given  by: 


uCx)  = - R 

It  should  be  noted  that  the  slope  du/dx  is  the  same  in  both  Regions  1 and 
5,  since  the  steady-state  mass  fluxes  at  the  lateral  boundaries  are  the 
same  in  both  regions. 

In  Regions  2 and  4 the  steady  state  source  terms  vanish,  and  the 
gas  velocity  becomes  a constant  given  by: 


u(x)  = + u^(  1 + B/R  ,1 


(55) 


where  the  upper  sign  refers  to  Region  2 and  the  lower  sign  refers  to  Region 
4.  The  magnitude  of  u(x)  was  obtained  by  substituting  x = 3/2  into  Equation 
(48).  The  particle  velocity  in  Region  2 is  then  given  by  (Appendix  D) : 


Up(x)  = u ( B/2  ) 1 


^-o(x-B/2) 

1+Tl 


(56) 


where  a = K/u(S/2K  For  x = B/2  the  exponential  term  is  unity  and  Equation 

(56)  yields  u (B/2)  = u(3/2)/(l  + ' ) which  is  the  velocity  of  the  particles 
P 

leaving  Region  1.  The  exponential  term  decays  rapidly  with  small  increases 
in  x to  yield 


u (x)  = u(  B /2)  = u,  ( 1 + B/R)  (57) 

P ° 

which  is  the  same  as  the  gas  velocity  in  Region  2.  The  particle  density 
quickly  relaxes  to  the  constant  value  given  by 

p (x)  - C (58) 

p m 

Again,  the  steady-state  particle  properties  in  Region  4 are  obtained  from 
the  Region  2 properties  by  symmetry. 

Center  Vent.  The  effect  of  a subsonic  exhaust  vent  at  the  center 

5 

of  the  T-burner  is  modeled  following  the  methodology  of  Levine  and  Culick. 
Here,  two  effects  are  considered:  (1)  the  mean  flow/acoustics  interaction, 

which  vanishes  in  the  absence  of  the  mean  flow,  and  (2)  the  acoustic  radiation 


32 


through  the  vent,  which  remains  finite  in  the  absence  of  mean  flow.  For 
odd  modes,  the  radiation  effect  is  weak  because  the  vent  is  located  at 
a pressure  node,  but  the  mean  flow/acoustics  effect  is  strong  due  to  the 
maximum  acoustic  velocity  at  this  location.  The  opposite  is  true  for  even 
modes;  the  mean  flow/acoustics  effect  is  weak  (velocity  node)  while  the 
radiation  effect  is  strong  (pressure  anti-node). 

In  considering  the  mean  flow/acoustics  vent  effect,  Levine  and 
Culick^  argued  that  the  net  effect  is  nearly  zero.  This  result  is  obtained 
because  the  vent  gain  effect  predicted  by  formal  application  of  linear 
theory  is  cancelled  by  a vent  loss  due  to  viscous  effects  and  the  axial 
force  exerted  by  the  walls  of  the  vent  on  the  flow.  Due  to  the  possibility 
of  a nonlinear  mean  flow/acoustics  vent  effect,  however,  Levine  and  Culick^ 
retained  the  ability  to  vary  the  mean  flow/acoustics  vent  effect  in  the 
nonlinear  differential  equations.  This  has  also  been  done  in  the  present 
approximate  analysis. 

In  the  vent  region  (Region  3 for  (1  - )/2  < x < (1  + 0^/2)), 

the  conservation  equations  (i.e..  Equations  (1)  - (5)),  are  modified  as 
follows.  The  source  terms  2m  /R  and  2tn  /R  are  replaced  by  id  and  <d  , respec- 
tively  which  are  defined  by; 


X 

P 


(60) 


where  A is  the  ratio  of  vent  area  to  chamber  cross-sectional  area,  and 

V 

u is  the  dimensionless  velocity  of  the  flow  out  the  vent.  The  minus  sign 
n 

indicates  that  x and  are  negative  for  flow  out  of  the  burner  (u^  positive), 
and  it  has  been  assumed  that  both  gas  and  particles  exit  with  the  same 
velocity  u . In  the  gas  and  particle  momentum  equations  (i.e.,  Equations 
(2)  and  (3))  the  axial  velocities  in  the  source  terms  arising  from  mean 
flow/acoustic  interactions  are  multiplied  by  the  factor  1 - . Setting 

=0  yields  vent  gain  which  corresponds  to  the  combustion  products  losing 
their  axial  momentum  as  they  exit  through  the  vent.  For  = 1 the  com- 
bustion products  retain  their  axial  momentum  and  the  net  vent  effect  is  zero, 
and  = 2 gives  a vent  loss.  Using  the  order  of  magnitude  analysis  given  in 
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Appendix  A,  the  modified  conservation  equations  are  combined  to  yield  the 
following  nonlinear  equations,  valid  in  the  vent  region; 


$ - i = 2ui  + (y+1-v  ) - — ^ 

XX  tt  xt  f,  dx  t 


+ + (Y-1) 


(%^t  ^ - a -V  (62) 


p'  = -Y  cp^  + ^ dj  ^ 


where 


uo  = - 


The  acoustic  radiation  vent  effect  is  described  by  the  term  u)^h^ 
appearing  in  Equation  (61).  This  introduces  the  additional’ unknown  quantity, 
u^  into  the  equations.  Assuming  plug  flow  in  the  pipe  connecting  the  burner 
with  the  surge  tank,  the  following  transient  vent  equation  is  obtained  from 
the  unsteady,  incompressible,  momentum  equation^: 
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'o  eff  dt 


Po  F = P 


In  this  equation,  effective  plug  length,  p ^ is  the  steady  state 

density,  is  the  surge  tank  pressure  (assumed  constant),  and  P is  the 
average  pressure  over  the  vent;  that  is. 


p(x',  t)dx' 


To  obtain  the  corresponding  equation  for  the  perturbation  u^  = u^  + u^ 

IS  introduced  into  Equation  (65)  and  the  steady  state  terras  are  subtracted 
out  to  give: 


eff  d 


^ + u u'  + 7(u')^  = f p'(x',  t)  dx' 

t n n 2 n 


where  p ^ = 1 has  been  used.  Neglecting  the  nonlinear  term  as  higher  order 
and  replacing  p'  by  its  first  order  equivalent  - gives 


du  ^ 

T n . - / 

L + u u 

eff  dt  n n 


^ ^^(x'.t)dx' 


In  the  vent  region  (Region  3)  the  steady-state  quantities  are 
also  obtained  by  integrating  the  steady-state  continuity  and  momentum 
equations,  where  the  steady-state  mass  source,  UD  , is  given  by: 
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(1  + 3/R) 


(69) 


U)  = 


2% 


V 


The  steady-state  gas  velocity  is  then  given  by: 

Uj^(l+3/R)  (l-2x) 
u(x)  = 


(70) 


for  (1  9y)/2  < X < (1  +3^)/2.  This  expression  shows  that  u(x)  decreases 

linearly  from  Uj^(l  + g/R)  at  the  left  boundary  of  Region  3 to  -u^(l  + 3 /R) 
at  the  right  boundary.  The  steady-state  gas  velocity  vanishes  at  the  center 
of  the  burner  (x  = 1/2).  The  particle  velocity  and  density  in  the  vent 
region  are  described  by  the  following  expressions  (Appendix  D); 


1 - 2x 


9| 


(71) 


Pp(x)  = 


1 - 


I.. 


l-2x 


q-1 


(72) 


where 


q = 2/Ti^  + (V^  - 1) 


and 


\ = 


+3/R) 


These  expressions  are  valid  for  (1  -3^)/2  < x < 1/2,  while  the  corres- 
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ponding  values  for  1/2  £ x < (1  +3^)/2  are  obtained  from  symmetry.  As 
in  Region  1 the  exponent  q is  large  for  particle  sizes  of  a few  microns; 
thus,  except  in  thin  zones  adjacent  to  the  boundaries  of  Region  3,  u^Cx) 
and  Pp(x)  are  well  approximated  by 


Up(x) 


(73) 


Pp(’^) 


(74) 


which  are  similar  to  Equations  (52)  and  (53)  for  Region  1.  Steady-state 
values  for  a typical  T-burner  configuration  with  cup  grains  are  shown  in 
Figure  4. 

An  additional  steady-state  value  is  needed  in  the  vent  region; 
that  is,  the  steady-state  velocity  of  the  flow  out  the  vent,  u^.  This 
quantity  appears  in  Equation  (68)  and  is  given  by: 


u = 
n 


2^  (1+e/R) 

A 


v 


(75) 


Amplitude  Equations  for  T-Burners.  To  summarize  the  solution  proce- 
dures for  T-burners,  the  application  of  the  Galerkin  method  in  the  five 
regions  previously  discussed  will  now  be  outlined.  The  series  expansions 
for  4 and  $ given  by  Equations  (19)  with  bj  = j'n  are  substituted  into 
the  appropriate  differential  equations  and  boundary  conditions  to  form 
the  residuals  E'^  T,  ^p}  > Regions  1, 

2,  4,  and  5,  Equations  (13)  and  (14)  are  used  to  obtain  E and  E^,  respec- 
tively, while  Equations  (61)  and  (62)  are  used  to  obtain  these  residuals 
in  Region  3.  The  Galerkin  orthogonality  conditions  given  by  Equations 
(46)  are  then  applied  to  obtain  a system  of  second-order  ordinary  differen- 
tial equations  for  the  mode-amplitudes  Aj(t).  The  spatial  integrations 
indicated  in  Equations  (46)  are  evaluated  piecewise  in  each  of  the  five 
regions,  using  the  appropriate  residuals  in  each  region. 

The  resulting  system  of  differential  equations  is  similar  to  that 


A 
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Dimensionless  Axial  Coordinate,  x 
Figure  4.  T-Burner  Steady-State  Properties 


obtained  for  motors  (i.e.,  Equations  (22)  and  (23)),  however,  some  differ- 
ences should  be  noted.  Due  to  the  acoustic  radiation  term  (i.e.,uu^h  ) 

c 

appearing  in  Equation  (61)  for  the  vent  region,  additional  terms  of  the 
form  C^(j)u'  appear  in  the  T-burner  mode-amplitude  equations.  Since  u^ 
is  an  additional  unknown  quantity,  another  equation  is  needed  in  order 
to  obtain  solutions.  The  required  equation  is  supplied  by  substituting 
the  series  expansion  for  § into  the  unsteady  flow  vent  equation  (i.e.. 
Equation  (68))  and  performing  the  indicated  integration  over  the  length 
of  the  vent.  The  resulting  equation  is  of  the  form; 

du^  dA  5 


where  u is  given  by  Equation  (75).  The  coefficients  of  the  linear  terms 
in  the  T-burner  equations  also  differ  from  the  corresponding  coefficients 
(i.e.,  Cj^  through  C^)  in  the  motor  equations,  due  primarily  to  the  differ- 
ences in  the  steady-state  properties  and  the  presence  of  burning  propellant 
disks  at  the  ends  of  the  burner.  Expressions  for  the  linear  coefficients 
in  the  T-burner  equations  are  given  in  Appendix  C.  The  coefficients  of 
the  nonlinear  terras  in  the  T-burner  equations  are  the  same  as  the  corres- 
ponding coefficients  in  the  motor  equations. 

2 . 7 Nonlinear  Combustion  Driving 

In  the  analysis  developed  in  the  preceding  sections,  the  response 

of  the  burning  rate  to  pressure  oscillations  was  described  by  a linear 

response  function  (see  Equations  (9)  through  (11)).  Based  on  intuition 

2 3 5 

and  studies  conducted  at  Georgia  Tech  ’ and  elsewhere  , it  is  expected 
that  as  the  amplitude  of  the  oscillation  grows,  the  true  combustion  re- 
sponse departs  fron  the  linear  value  predicted  by  Equation  (9).  Thus,  a 
nonlinear  expression  relating  the  perturbations  in  burning  rate  to  the 
gas  flame  pressure  fluctuations  is  needed,  if  accurate  stability  predic- 
tions are  to  be  obtained. 

A nonlinear  pressure  coupled  response  is  presently  included  in 
the  "exact"  analysis  developed  by  Kooker  and  Zinn  (in  fact  one  of  the  mod- 
ifiaacions  of  the  "exact"  theory  required  for  the  study  of  gasdynamic 
non-linearities  was  to  replace  this  nonlinear  response  with  a linear  one). 
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This  nonlinear  response  is  obtained  by  solving  the  governing  equations 

describing  the  thermal  wave  in  the  solid,  the  surface  pyrolysis  reactions, 

and  the  gas  phase  reactions  simultaneously  with  the  conservation  equations 

desc'ibing  the  unsteady  flow  in  the  chamber  and  nozzle.  In  order  to  use 

the  same  approach  in  the  approximate  analysis,  it  would  be  necessary  to 

apply  the  Galerkin  Method  to  obtain  the  behavior  of  the  temperature  in 

the  solid  propellant  and  the  flame  zone.  In  principle,  this  solution  for 

the  combustion  response  would  have  to  be  obtained  simultaneously  with  the 

solutions  of  the  combustor  conservation  equations,  as  was  done  in  the 
2 3 

"exact"  analyses  ’ Attempts  to  use  the  Galerkin  Method  to  obtain  the 
nonlinear  solid  propellant  combustion  response  were  made  during  the  develop- 
ment of  the  Kooker-Zinn  model;  however,  these  studies  have  shown  that  the 
approximate  technique  was  unsuitable  for  properly  predicting  the  nonlinear 
propellant  response.  Based  upon  this  experience,  it  is  expected  that  the 
same  difficulties  would  also  be  encountered  when  attempting  to  incorporate 
the  nonlinear  combustion  response  into  the  present  approximate  analysis. 
Thus,  in  the  absence  of  any  other  nonlinear  models  for  pressure  coupling, 
a heuristic  combustion  model  has  been  used  in  the  studies  to  assess  the 
importance  of  modelling  the  nonlinear  combustion  response. 

A number  of  possible  expressions  for  describing  the  nonlinear  burn- 
ing rate  could  be  obtained  by  adding  various  nonlinear  terms  to  the  linear 

2 

burning  response,  such  as  terras  proportional  to  | p'  |,  (p')  , | p'l  (p‘), 

3 

and  (p')  . Since  the  primary  concern  here  is  determining  the  effect  of 

nonlinearities  in  the  combustion  response  upon  the  nonlinear  stability 

characteristics  of  a motor,  rather  than  determining  the  correct  form  of 

the  nonlinear  combustion  response,  the  simplest  of  these  heuristic  models 

was  used.  Thus  a nonlinear  combustion  response  was  incorporated  into  the 

approximate  analysis  by  replacing  the  linear  response  functions  R appear- 

m 

ing  in  Equations  (22)  with  nonlinear  response  functions  given  by; 

m 


' 

dA 

■■ 

1 + b 

m 

m 

dt 

_ 

(77) 


In  Equation  (77)  is  the  linear  response  factor  for  the  m*"^  mode,  b^ 

is  a complex  constant,  and  | dA^/dt  | is  the  magnitude  of  the  complex  arapli- 
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tude  function  given  by 


where  the  superscripts  r and  i denote  real  and  imaginary  parts  of  A^,  re- 
spectively. Since  the  pressure  amplitude  is  proportional  to  [dA/dt] , the 
effective  nonlinear  response  factor,  is  a complex  number  whose  magni- 

tude and  phase  are  linear  functions  of  the  amplitude  of  the  oscillation. 

This  expression  reduces  to  the  linear  burning  rate  expression  in  the  limit 
of  small  amplitude  oscillations.  At  larger  amplitudes  the  magnitude  of 
the  response  factor  will  be  increased  or  decreased  from  its  linear  value 
according  to  the  value  of  b;  also,  the  phase  will  be  shifted  if  the  imagin- 
ary part  of  b is  nonzero.  Calculated  solutions  obtained  using  this  heuris- 
tic nonlinear  combustion  model  are  presented  in  section  4.3,  where  they 
are  compared  with  approximate  solutions  obtained  with  the  linear  combustion 

response  (b  = 0)  as  well  as  with  "exact"  solutions  obtained  using  the 
m 

Kooker-Zinn  nonlinear  combustion  model. 

2 . 8 Nonlinear  Particle  Damping 

In  the  preceding  sections  it  was  assumed  that  the  viscous  interaction 
between  the  particles  and  gas  was  described  by  the  Stokes  Drag  Law  (Equations 
(7)  and  (8))  for  which  the  drag  force  is  proportional  to  the  relative  velo- 
city between  the  particles  and  gas.  This  linear  drag  law  is  valid  for 
small  Reynolds  numbers  up  to  about  unity.  For  many  situations  encountered 
in  solid  rocket  motors,  however,  the  particle  sizes  and  frt.  , ncies  are 
sufficiently  large  that  higher  Reynolds  numbers  can  occur.  In  these  situa- 
tions the  Stokes  law  is  no  longer  applicable  and  a nonlinear  drag  law  must 

be  used.  The  nonlinear  particle  drag  effects  are  described  by  a higher 

21 

order  correction  to  the  Stokes  law  given  by 
21 

Rudinger,  G.,  "Effective  Drag  Coefficients  for  Gas-Particle  Flow 
in  Shock  Tubes,"  Project  Squid  Technical  Report  CAL-97-PU,  March  1969. 


41 


(79) 


Fp  = Kpp(u  - Up) 


1 + 


where  Re  is  the  Reynolds  number  given  by 


Re 


D*  c*  cj* 
o 


P 1 u - u 1 


(80) 


The  straightforward,  rigorous  method  of  introducing  the  nonlinear 
particle  drag  law  (Equation  (79))  into  the  approximate  analysis  is  to  sub- 
stitute Equation  (79)  for  Fp  where  it  appears  in  the  gas  and  particle 
momentum  equations  (i.e.,  Equations  (3)  and  (4))  and  the  energy  equation 
(i.e..  Equation  (5))  and  to  follow  the  same  order  of  magnitude  analysis 
outlined  in  Section  2.2  and  Appendix  A.  This  procedure  leads  to  insurmount- 
able difficulties  due  to  the  highly  nonlinear  form  of  the  drag  correction 
term  (i.e.,  the  absolute  value  of  the  relative  velocity  raised  to  the  2/3 
power).  The  first  difficulty  arises  when  attempting  to  separate  the  steady- 

state  component  from  the  equations.  The  relative  velocity  term  can  be 

• " 2/3 1 • • 1 2/3 

written  as  1 u - u | 1 1 + (u'  - u'  )/(u  - u )|  , but  a binomial  expan- 

F.  . .F  P'-_ 

sion  of  this  term  is  valid  only  if  (u'  - u'  )/(u  - u )<<1.  Since  the 

P P 

relative  velocity  perturbation  u'  - Up  may  often  be  of  the  same  order  of 
magnitude  as  the  steady  state  relative  velocity  u - Up,  the  binomial  ex- 
pansion can  not  be  used  and  the  steady  state  and  perturbation  quantities 
can  not  be  separated.  Even  if  this  difficulty  did  not  exist,  other  problems 
arise  when  attempting  to  combine  the  conservation  equations  to  obtain  the 
gas  and  particle  potential  equations  by  the  procedure  followed  in  Appendix 
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A.  Finally,  the  absolute  value  in  the  drag  correction  terms  make  it  imposs- 
ible to  separate  the  time  and  space  variables  when  applying  the  Galerkin 
method,  a crucial  step  needed  in  order  to  evaluate  the  spatial  integrations 
involved. 

Because  of  these  difficulties  in  obtaining  equations  for  the  gas 
and  particle  potentials  with  given  by  Equation  (79),  a heuristic  approach 
is  followed.  It  is  tacitly  assumed  that  the  effect  of  particle  drag  non- 
linearities  can  be  described  by  the  previously  derived  mode-amplitide  equa- 
tions (i.e..  Equations  (22)  and  (23))  in  which  the  linear  drag  constant 
K is  replaced  by  a nonlinear  drag  coefficient  K which  is  amplitude  depen- 
dent. It  is  postulated  that  K has  the  following  form: 


1 + K + CK„, 
ss  NL 


A - A 


2/3 


(81) 


where  K is  the  original  linear  drag  constant  given  by  Equation  (8), 
is  a steady-state  contribution  to  the  nonlinear  drag  correction,  is 

given  by 


K 


NL 


1 t 

M I 


(82) 


and  C is  an  adjustable  parameter.  The  steady-state  contribution  is  ob- 
tained by  averaging  the  nonlinear  drag  term  corresponding  to  the  steady- 
state  relative  velocity,  u(x)  - u^(x),  over  the  length  of  the  motor.  Since 
the  steady-state  relative  velocity  is  proportional  to  x (Equations  (40) 
and  (41)),  the  value  of  is  given  by : 
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ss 


i i 
6 1 


c*  o* 
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2/3 
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,2/3 
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2/2. 

X dx 


(83) 


Evaluating  the  integral  and  using  Equations  (41)  and  (82)  give 
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K = I K„, 

ss  5 NL 


1 + 8u  /K  - 1 


1 +Jl  + 8u^/K 


The  unsteady  contribution  to  the  drag  correction  is  assumed  to  be  propor- 

tional  to  the  quantity  |A  - A^|  where  A and  A^  are  the  complex  gas  and 

particle  mode -amplitude  functions  for  the  fundamental  mode  (the  effect 

of  the  higher  modes  upon  this  term  has  been  neglected).  The  constant  of 

2/3 

proportionality  for  this  term  contains  the  O dependence  expected  from 
Equation  (79),  and  the  adjustable  parameter  accounts  for  the  axial  varia- 
tion of  the  relative  velocity  perturbation  u'-  u\ 

The  heuristic  nonlinear  particle  drag  model  described  above  has 
been  incorporated  into  the  approximate  computer  program.  Calculated  solu- 
tions obtained  with  this  model  are  presented  in  Section  4.4  where  they 
are  compared  with  the  available  "exact"  solutions.  On  the  basis  of  these 
comparisons  the  validity  of  the  heuristic  nonlinear  particle  drag  model 
is  assessed. 


'•  ' r- 
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3.  EXACT  ANALYSIS 


2 3 

The  "exact"  analysis  developed  by  Kooker  and  Zinn,  ’ which  was 
used  in  the  present  investigation  is  based  on  a finite  difference  solution 
of  the  quasi-one-dimensional  conservation  equations  for  a solid  rocket 
motor  and  nozzle  combination.  These  equations  are  given  as  follows: 


Continuity  (gas) 


3p  , du  ^ Sp 

9t  3x  9x 


Continuity  (particles) 


- + p r + u r — 

dt  p ox  p ox 


p dlnA 

R ■ Pp^p  ~ 


Momentum  (gas) 


Su  , 9u  ^ 1 dp 

ot  Yp 


= - |Kpp(u  - Up)  + f mgU  }/p 


Momentum  (particles) 


Energy  (gas) 


r — + u r — 

ot  p dx 


= K(u  - u ) - f 


= 

dt  dx 


(tw^) 


where 


Umu^  + m 1 + Kp  (u  - u )^1 

2 \ g P P / J ^P  P I 


S = Inp  - Y Itip 


The  quantities  m and  ra  respectively  represent  the  mass  flow  rates  of 
g P 

gas  and  particles  per  unit  surface  area  leaving  the  burning  propellant  surface. 
Here  all  mass  addition  due  to  the  burning  propellant  occurs  at  the  combustor 
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boundaries  and  the  volume  sources  due  to  evaporation  or  combustion  of  parti- 
cles are  neglected.  The  source  terms  in  the  momentum  equations  result 
from  the  application  of  Stokes'  Law  to  describe  the  particle  drag  and  the 
force  required  to  accelerate  the  incoming  gas  and  particles  to  their  local 
velocities.  The  source  terms  in  the  energy  equation  arise  from  the  local 
addition  of  hot  combustion  products  (with  enthalpy  h^ ) entering  the  control 
volume  from  the  flame  zone  and  particle  drag  within  the  control  volume. 
Including  the  area  variation  term  in  the  continuity  equation  allows  a quasi- 
one-dimensional  solution  for  the  unsteady  flow  in  the  nozzle  as  well  as 
in  the  combustion  chamber.  A more  complete  description  of  the  combustor 
model  and  method  of  numerical  integration  of  the  equations  is  given  in 
Chapter  II  of  Reference  3. 

To  complete  the  analysis,  an  unsteady  combustion  model  is  needed 
to  determine  the  time-dependent  quantities  m^  and  m^  which  appear  as  mass 
and  energy  sources  at  the  boundary  of  the  combustor.  This  is  done  in  the 
Kooker-Zinn  model  by  solving  the  transient  burning  rate  response  to  the 
imposed  chamber  pressure  oscillations  simultaneously  with  the  chamber  conser- 
vation equations.  The  combustion  model  used  at  each  burning  station  along 

the  propellant  is  based  on  the  following  four  simplifying  assumptions  which 

12 

have  been  employed  in  nearly  all  combustion  instability  studies  to  date: 

(1)  the  unburned  propellant  is  homogeneous  and  one-dimensional,  (2)  the 
conversion  of  the  solid  phase  to  gas  is  represented  by  a simple  Arrhenius 
type  pyrolysis  law,  (3)  condensed  phase  reactions  are  neglected,  (4)  the 
behavior  of  the  gas-phase  flame  zone  is  quasi-steady.  Under  these  assump- 
tions the  problem  separates  into  three  distinct  regions  which  are  analyzed 
individually:  (1)  the  nonreacting  solid  phase  region,  (2)  the  solid-gas 

interface  where  pyrolysis  occurs,  and  (3)  the  gas-phase  flame  zone.  The 
conservation  equations  describing  the  behavior  of  this  model  and  their 
numerical  solutions  are  discussed  in  Chapter  III  of  Reference  3. 

For  the  case  of  a constant  area  chamber,  the  above  equations  are 
identical  to  Equations  (1)  through  (6)  upon  which  the  approximate  analysis 
is  based.  Thus  both  the  "exact"  and  approximate  analyses  are  attempting 
to  solve  the  same  set  of  equations.  However,  the  combustion  response  model 
and  the  treatment  of  the  nozzle  used  in  the  original  Kooker-Zinn  model 
are  different  from  those  used  in  the  approximate  analysis.  In  order  to 
facilitate  comparisons  with  the  approximate  analysis,  the  original  Kooker- 
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Zinn  model  was  modified  in  the  present  study.  A brief  description  of  the 
required  modifications  is  given  below. 


3 . 1 Modification  for  Quasi-Steady  Nozzle 

In  the  original  Kooker-Zinn  model  the  nonlinear  governing  equations 
were  solved  numerically  throughout  the  entire  system  consisting  of  both 
the  combustion  chamber  and  the  nozzle.  The  downstream  boundary  condition 
was  imposed  at  an  arbitrary  location  downstream  of  the  nozzle  throat  where 
the  flow  is  supersonic.  The  effect  of  the  nozzle  obtained  from  such  an 
analysis  is  inherently  nonlinear  and  would  be  expected  to  differ  significantly 
from  the  nozzle  effect  obtained  from  the  linear  nozzle  admittance  boundary 
condition  used  in  the  approximate  analysis.  Therefore  the  "exact"  analysis 
was  modified  by  eliminating  the  nozzle  form  the  flowfield  and  replacing 
it  with  a quasi-steady  linear  nozzle  admittance  condition  to  be  satisfied 
at  the  nozzle-end  of  the  chamber.  The  method  of  characteristics  is  used 
to  obtain  the  boundary  values  at  the  nozzle-end  in  a manner  similar  to 
that  used  at  the  head-end.  This  modification  allows  the  nozzle  effect 
to  be  the  same  in  both  the  approximate  and  "exact"  analyses. 

At  the  nozzle-end  of  the  chamber,  the  velocity  and  pressure  must 
satisfy  the  following  relation: 
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[ where  Y is  the  nozzle  admittance  and  u^  and  p^  are  the  corresponding  steady- 

1 state  quantities.  For  a quasi-steady  nozzle  (neglecting  the  effect  of 

, particles)  the  nozzle  admittance  Y is  given  by: 
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Y -1  e 

" Pe 


(92) 


These  relations  have  been  incorporated  into  the  method  of  characteristics 
solution  for  the  dependent  variables  at  the  nozzle-end  boundary.  The  nozzle 
calculations  are  performed  by  a new  subroutine  (NOZMOC)  which  was  added 
to  the  existing  "exact"  computer  program. 

A test  case  was  run  at  this  point  to  check  out  the  quasi-steady 
nozzle  modification.  Pressure  waveforms  calculated  with  a conventional 
(i.e.,  nonlinear)  nozzle  and  the  quasi-steady  nozzle  are  compared  in  Figure 
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5.  Here  it  is  seen  that  substituting  the  quasi-steady  nozzle  for  the  conven- 
tional nozzle  results  in  an  increase  in  the  frequency  of  the  oscillation 
and  a decrease  in  the  decay  rate.  The  increased  frequency  results  from 
the  decrease  in  length  over  which  the  wave  travels,  and  the  decrease  in 
damping  results  because  the  quasi-steady  nozzle  is  more  like  a hard-wall 
termination  than  the  conventional  nozzle. 

3 . 2 Linearized  Burning  Rate  Model 

The  other  modification  of  the  "exact"  analysis  concerns  the  transient 
burning  rate  model.  In  the  original  "exact"  analysis  solutions  were  obtained 
for  the  thermal  profiles  in  the  solid  propellant,  which  yield  the  surface 
temperature  and  regression  rate  at  each  increment  of  time.  These  profiles 
were  coupled  to  the  pressure  oscillations  in  the  chamber  through  the  reaction 
kinetics  in  the  flame  zone  and  at  the  solid  surface.  These  processes  yielded 
a nonlinear  response  of  the  burning  rate  to  the  pressure  oscillations. 

Since  the  above  model  cannot  be  readily  incorporated  into  the 
approximate  analysis,  the  approximate  solutions  were  obtained  using  a linear 
combustion  response  function  Equation  (9).  In  order  to  compare  the  results 
obtained  with  the  approximate  and  "exact"  theories,  a linear  combustion 
response  model  was  incorporated  into  the  "exact"  analysis.  Since  the  fre- 
quency and  waveshape  of  the  pressure  oscillation  are  not  known  a priori, 
the  two-parameter  (A,B)  form  of  the  linearized  response  function  (see  Equation 
(10))  could  not  be  introduced  directly  into  the  "exact"  program.  Instead, 
a perturbation  analysis  was  performed  on  the  nonlinear  equations  governing 
the  response  of  the  burning  solid  propellant  to  pressure  disturbances, 
and  the  resulting  higher  order  (i.e.,  nonlinear)  terms  were  neglected. 

The  resulting  linear  equations  and  boundary  conditions  are  given  below: 

Energy  Equation  in  Solid  Propell; 

T^ ' + r T 

Surface  Regression  Law 
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Pressure  (atmospheres) 


Figure  5.  Comparison  of  Pressure  Waveforms  for  Conventional 
and  Quasi-Steady  Nozzles 
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Boundary  Conditions 


x'(-y^.t)  = 0 


t;  (o.t)  . -Zit'  - f (cVc*  - l)  I' 
f [ p r _ 


In  the  above  equations  the  barred  quantities  are  steady  state  values,  while 

the  primed  quantities  are  perturbations.  The  quantities  p,  r,  and  T(y) 

respectively  denote  pressure  at  the  propellant  surface,  surface  regression 

rate,  and  temperature  in  the  solid  propellant  at  depth  y.  The  quantities 

E^,  Z^,  and  are  constants  which  depend  on  the  propellant  and  gas  properties. 

A computer  subroutine  (LINTHW)  was  developed  to  solve  the  linearized 

equations  given  above.  This  subroutine,  which  is  based  on  the  method  of 

invariant  imbedding  used  in  the  original  nonlinear  program  developed  by 

Kooker,  calculates  the  linearized  response  of  the  surface  regression  rate 

r (and  hence  the  mass  flux  m ) to  small  amplitude  pressure  disturbances 

g 

of  arbitrary  waveform.  For  sinusoidal  pressure  disturbances  the  linearized 

response  computed  by  LINTHW  is  equivalent  to  the  two-parameter  (A,B)  model. 

In  order  to  compare  "exact"  and  approximate  solutions,  values  of 

A,  B,  and  n corresponding  to  a given  set  of  propellant  properties  are  needed. 

These  are  obtained  by  an  analytical  solution  of  the  above  linearized  combustion 

response  model  for  assumed  sinusoidal  pressure  disturbances.  Explicit 

relations  for  A and  B in  terms  of  the  solid  propellant  activation  energy 

E*,  the  steady  state  surface  temperature  T*,  and  the  surface  heat  of  reaction 
s s 

Q*  are  easily  derived  and  are  given  by: 
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In  the  linear  analysis,  the  value  of  n cannot  be  obtained  explicitly 
in  terms  of  propellant  properties  but  must  be  obtained  from  burning  rate 
calculations  at  low  frequencies.  An  auxiliary  program  which  calculates 
the  burning  rate  corresponding  to  a sinusoidal  pressure  input  has  been 
developed  to  determine  the  complex  response  factor  R by  comparing  the 
computed  amplitude  and  phase  of  the  burning  rate  perturbation  with  that 
of  the  pressure  perturbation.  This  program  is  then  used  to  determine  the 
value  of  n,  since  the  computed  values  of  ft  approach  n as  frequency 
approaches  zero.  After  n is  determined,  values  of  ft  are  obtained  for  several 
different  frequencies  using  both  the  linearized  and  nonlinear  burning  rate 
models.  These  values  are  compared  with  the  response  curve  obtained  analy- 
tically from  Che  two-parameter  linear  response  function  described  by  Equa- 
tions (10)  and  (11).  An  example  of  such  a comparison  is  shown  in  Figure 
6 where  the  close  agreement  obtained  is  an  indication  of  the  validity  of 
both  the  linearized  and  nonlinear  combustion  models  for  small  amplitude 
oscillations . 

3 . 3 Nonlinear  Particle  Damping 

In  the  original  Kooker-Zinn  "exact"  model  the  particle  drag  was 
described  by  the  Stokes  Drag  Law,  which  is  a linear  expression  valid  for 
Reynolds  numbers  of  order  unity.  In  order  to  include  nonlinear  particle 

drag  effects  which  occur  at  higher  Reynolds  numbers  (larger  particle  sizes 

21 

and  higher  frequencies),  the  higher  order  correction  to  the  Stokes  Law 
was  incorporated  into  the  "exact"  analysis.  Thus  the  drag  term  in  the 
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momentum  equation  is  given  by  Equations  (79)  and  (80).  A more  convenient 
form  for  numerical  calculations  is  given  by 


F = Kp  (u 
P P 


-v[ 


1 + V 2/3, 

1 + |u-u 


where  is  given  by  Equation  (82).  These  expressions  are  readily  incor- 

porated into  the  appropriate  conservation  equations  (i.e..  Equations  (87) 
through  (89))  by  replacing  the  drag  constant  K with  the  variable  drag 
coefficient  K defined  by: 


K = K 1 + 


The  above  nonlinear  drag  law  has  been  incorporated  into  the  Kooker-Zinn 
"exact"  computer  program;  the  linear  drag  law  option  is  available  by  speci- 
fying = 0. 


4.  RESULTS  AND  DISCUSSION 


In  this  section,  numerical  calculations  of  nonlinear  axial  mode 
instabilities  in  solid  rocket  motors  and  T-burners  are  presented  and  dis- 
cussed. Most  of  the  approximate  solutions  presented  herein  were  calculated 
using  the  Galerkin  method  without  averaging;  that  is,  by  obtaining  numerical 
solutions  of  the  system  of  second  order  differential  equations  given  by 
Equations  (22)  and  (23).  These  approximate  solutions  are  compared  with 
"exact"  solutions  obtained  by  the  Kooker-Zinn  model  in  many  cases.  The 
approximate  solutions  are  also  compared  with  Culick's  approximate  solutions 
(Method  of  Averaging)  and  Levine's  "exact"  solutions  when  available. 

Most  of  the  approximate  solutions  presented  in  this  section  were 
obtained  in  a parametric  study  to  assess  the  importance  of  gasdynamic  non- 
linearities  (mode-coupling)  on  the  stability  of  solid  rocket  motors,  a major 
objective  of  this  program.  In  this  study,  the  only  nonlinear  process  consid- 
ered was  gasdynamical  mode-coupling,  while  other  processes  such  as  combus- 
tion driving  and  particle  damping  were  described  by  linear  models.  The 
following  parameters  were  considered:  (1)  the  number  of  modes  used  to 

describe  the  total  wave,  (2)  the  magnitude  and  harmonic  content  of  the 
initial  disturbance,  (3)  the  characteristics  of  the  propellant  response 
function,  (4)  the  oscillation  frequency,  and  (5)  the  size  and  concentration 
of  the  aluminum  oxide  particles.  The  approximate  and  exact  analyses  are 
comsared  on  the  basis  of  the  predicted  growth  or  decay  rates  for  the  tran- 
sient solutions  and  the  final  limiting  amplitude  and  waveform.  In  addition, 
the  transient  behavior  and  relative  amplitudes  of  the  individual  modes  are 
determined . 

In  addition  to  the  above  parametric  study  approximate  solutions 
are  presented  for  (1)  motors  with  a nonlinear  combustion  response  to  pressure 
oscillations,  (2)  motors  with  nonlinear  particle  damping,  and  (3)  T-burners. 
Also  solutions  obtained  using  the  Method  of  Averaging  (MOA)  are  compared 
with  the  Galerkin  method  solutions  for  motors  with  and  without  particles. 
Finally  approximate  solutions  are  compared  with  available  experimental  data 
for  both  motors  and  T-burners. 

4 . 1 Typical  Nonlinear  Solutions 

Betore  presenting  the  results  of  the  parametric  studies,  the  applica- 
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tion  of  the  approximate  nonlinear  analysis  will  be  illustrated  by  an  example. 

To  obtain  an  approximate  solution  for  a given  engine  configuration 
and  operating  condition,  the  following  must  be  specified;  (1)  the  number 
of  modes  N present  in  the  approximate  series  solution;  (2)  the  itructure 
of  the  initial  disturbance  (i.e.,  A.(0)  and  dA^(0)/dt  for  j = 1,  2,  ...  2N); 
(3)  the  combustor  and  propellant  grain  lengths;  (4)  the  nozzle  admittance 
Y for  each  of  the  modes  in  the  solution;  (5)  the  combustion  parameters  A, 

B,  and  n or  the  values  of  R for  each  mode;  (6)  the  concentration  and  size 
of  the  aluminum  oxide  particles;  (7)  the  mean  flow  Mach  number  at  the  nozzl-- 
entrance;  (8)  the  steady-state  temperature  in  the  combustor  and  (9)  the 
ratio  of  the  specific  heats  of  the  gaseous  phase.  For  the  example  consid- 
ered here  a five-mode  series  was  employed,  the  initial  disturbance  was  a 
pure  fundamental  mode  (IL)  oscillation  of  about  15%  pressure  amplitude, 
and  the  nozzle  admittances  for  each  of  the  modes  were  obtained  by  the  quasi- 
steady nozzle  relation  (i.e..  Equation  (43);.  The  specific  heat  ratio  was 
Y = 1.2  and  the  mean  flow  Mach  number  at  the  nozzle  entrance  was  M = 0.10. 
The  response  function  was  described  by  A = 7.54,  B = 0.686,  and  n = 0.81 
(Figure  7).  This  is  the  moderately  strong  response  considered  by  Kooker 
in  the  "exact"  analysis.  Here  the  frequency  fi  was  chosen  such  that  all 
of  the  modes  lay  on  the  descending  branch  of  the  response  curve;  thus,  the 
real  parts  of  the  responsj  function  for  the  higher  frequency  modes  were 
all  smaller  than  the  corresponding  value  for  the  fundamental  (IL)  mode. 

For  such  cases  the  series  expansion  in  terms  of  the  acoustic  modes  was  expec- 
ed  to  converge  fairly  rapidly,  and  the  mode-coupling  should  tend  to  limit 
the  amplitude.  Finally,  it  was  assumed  that  no  particles  were  present  in 
the  combustor . 

The  approximate  solutions  (Galerkin  method)  for  this  case  are  shown 
in  Figures  8 and  9.  Figure  8 shows  the  relative  magnitudes  and  transient 
development  of  the  individual  modes  (i.e.,  the  real  parts  of  the  functions 
A.(t))  following  the  introduction  of  the  pure  IL  mode  initial  disturbance. 

It  is  seen  that  the  amplitudes  of  the  higher  harmonics  are  much  smaller 
than  the  amplitude  of  the  fundamental  mode.  The  transient  development  of 
the  head-end  pressure  perturbation  shown  in  Figure  9 was  determined  from 
the  individual  modes  using  Equations  (19)  and  (24).  The  initially  sinu- 
soidal pressure  waveform  quickly  beccmes  distorted  due  to  the  excitation 
of  the  higher  harmonics. 
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Dimensionless  Frequency,  C 


Figure  7. 


Moderate  Response  Curve  with  Positions  of  First  Five 
Axial  Modes 


Real  Part  of  Mode -Amplitude  Function 
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Figure  9,  Head-End  Pressure  Waveform  for  Unstable  Motor 


It  should  be  noted  that  the  individual  mode-amplitude  functions 
in  Figure  8 are  the  real  parts  of  the  A^Ctl's  for  the  gas-phase  velocity 
potential;  thus,  they  are  the  harmonic  components  of  the  velocity  waveform 
rather  than  the  pressure  waveform.  From  Equations  (19a)  and  (24)  and  the 
nearly  sinusoidal  shape  of  the  Aj(t)'s,  it  is  easily  seen  that  the  amplitudes 
of  the  corresponding  harmonics  of  the  pressure  waveform  are  approximately 
obtained  by  multiplying  the  amplitudes  of  the  Aj(t)'s  by  yjoui  where  aj  j 
is  the  frequency  of  the  IL  mode.  This  relationship  is  valid  for  oscillations 
of  sufficiently  small  amplitude  that  the  second  order  terms  in  equation 
(24)  are  negligible.  For  larger  oscillation  amplitudes  the  individual 
pressure  harmonics  may  be  obtained  by  Fourier  analysis  of  the  waveform 
obtained  by  Equation  (24). 

4 . 2 Parametric  Studies  of  Mode-Coupling 

The  basic  set  of  motor  parameters  used  in  this  study  was  based  on 

3 

the  data  given  by  Levine  and  Culick  which  correspond  to  a small  laboratory 
pulse  motor.  This  motor  has  a grain  length  of  about  .597  m (1.958  ft)  and 
an  initial  bore  of  about  50  mm  (2  in).  The  physical  properties  of  the 
propellant  (ANB  3066),  and  its  gaseous  and  particulate  combustion  products, 
are  given  in  Table  1. 

The  data  along  with  the  endothermic  heat  release  at  the  propellant 

surface  (Qg)  ts  used  as  input  to  the  Kooker-Zinn  "exact"  analysis,  which 

then  computes  the  corresponding  values  of  A,  B,  and  n which  describe  the 

transient  propellant  response.  A reference  state  is  also  needed  in  order 

to  evaluate  some  of  the  constants  which  appear  in  the  burning  rate  model; 

in  this  case  a surface  regression  rate  of  1.17  cm/sec  (0.459  in/sec)  and 

a surface  temperature  of  800°K.  (1440*^R)  is  assumed  for  a pressure  of  10810 
2 

knt/m  (1568  psi).  For  most  cases  considered  here,  the  mean  chamber  pres- 
sure is  assumed  equal  to  the  above  reference  pressure.  For  this  case  the 
total  flux  of  burned  gases  leaving  the  tubular  propellant  grain  is  2.338 
kg/sec  which  corresponds  to  a Mach  number  at  the  nozzle  entrance  of  = 
0.0780.  For  this  combination  of  chamber  length  and  steady-state  temperature 
(Tj*),  the  chamber  sonic  speed  is  1278  m/sec  (4194  ft/sec)  giving  a funda- 
mental mode  frequency  of  1071  Hz  for  the  pure  gas  (without  particles) 

When  comparing  the  approximate  solutions  with  the  "exact"  solutions, 
the  values  of  A,  B,  and  n computed  from  the  above  data  are  used  as  input 
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Table  1.  ANB  3066  Propellant  Data 


Propellant 


densi ty : 

specific  heat: 

c* 

s 

thermal  conductivity: 

k* 

s 

activation  energy: 

E* 

s 

Gas-Phase  Flame 

specific  heat  ratio: 

y 

specific  heat 

c* 

p 

thermal  conductivity 

k* 

g 

activation  energy: 

^*f 

flame  temperature: 

Particles 

density:  p* 

m 

diameter:  O’ 

C 

m 


3 

= 1,766  gm/cra 

= 0.329  cal/gm-°K 

= 1,791x10  ^ cal/cm-sec- 

= 15,250  cal/mole 

1.23 

= 0.483  cal/gm-°K 

-4 

= 1.658x10  cal/cm-sec- 

*=  30,000  cal/mole 

= 3525°K  (6345°R) 

= 4,0  gm/cm 

2.5  n 


particle  loading: 


0.36 


to  the  approximate  analysis.  In  order  to  completely  define  the  values  of 
and  for  each  of  the  modes,  the  frequency  parameters  are 

also  needed  for  each  mode.  These  are  obtained  from  the  computed  frequency 
of  the  IL  mode  and  the  steady  state  burning  rate  obtained  from  the  "exact" 
analysis,  and  are  given  by: 

j = 1,2, ...N  (101) 

The  most  useful  parameter  to  vary  to  obtain  different  response  functions 

^ “it 

is  Q ; making  Q more  negative  (exothermic)  decreases  B and  yields  larger 
values  of  the  peak  response.  The  frequency  parameter  is  varied  by  changing 
the  chamber  length  or  by  changing  the  steady  state  burning  rate. 

The  results  of  the  parametric  studies  will  now  be  presented  in  the 
following  order:  (1)  effect  of  number  of  modes  in  series,  (2)  effect  of 
initial  disturbance,  (3)  effect  of  propellant  response  function,  (4)  effect 
of  frequency,  and  (5)  effect  of  particles. 

Effect  of  Number  of  Modes  Used  to  Represent  the  Solution.  Figures 
10  and  11  show  the  effects  of  the  number  of  modes  N used  to  represent  the 
approximate  solutions  for  a case  in  which  no  particles  are  present.  Solu- 
tions were  obtained  using  a basic  series  consisting  of  the  three  lowest 
frequency  axial  modes  (i.e.,  the  first,  second  and  third  longitudinal  modes). 
Higher  harmonics  were  then  added  one  at  a time  to  the  basic  series  and  the 
resulting  solutions  were  compared  on  the  basis  of  wave  shape  and  limiting 
amplitude.  This  comparison  was  used  to  determine  the  minimum  number  of 
modes  needed  to  adequately  describe  the  behavior  of  the  unstable  rocket 
motor  with  a minimum  expenditure  of  computer  time. 

Figure  10  shows  the  development  of  the  peak-to-peak  head-end  pressure 
amplitude  with  time  for  an  initial  7%  fundamental  mode  disturbance.  The 
combustion  parameters  are  A = 6.00,  B = 0.590,  n = 0.583  and  = 4.907. 

The  growth  of  the  pressure  oscillation  amplitude  is  shown  for  one,  three 
four,  five,  and  six  mode  series  expansions.  For  the  one-mode  series 
(N  = 1)  expansion,  there  is  no  nonlinear  mode-coupling  to  limit  the  ampli- 
tude, so  the  oscillations  grow  at  the  linear  rate  of  about  80  sec  . For 
4 < N < 6,  increasing  the  number  of  modes  decreases  the  predicted  limit- 
cycle  amplitude  as  a result  of  the  nonlinear  coupling  of  the  fundamental 
mode  with  increasingly  stable  higher  modes.  The  increased  stability  of 
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Peak- to- Peak  Head-End  Pressure  Amplitude 


these  higher  modes  results  from  the  decreased  combustion  response  at  higher 
frequencies  for  this  case  (Table  2).  As  the  number  of  modes  is  increased, 
the  approximate  solutions  also  approach  the  "exact"  solution  also  shown 
in  Figure  10,  but  it  is  not  clear  from  this  figure  that  the  approximate 
solutions  converge  to  the  "exact"  solution  in  the  limit  of  an  infinite 
number  of  modes.  Indeed,  it  appears  that  the  approximate  solutions  will 
converge  to  a larger  limiting  amplitude  than  the  "exact"  solutions,  probably 
due  to  the  order  of  magnitude  approximations  made  in  deriving  Equations 
(13)  and  (14)  upon  which  the  approximate  technique  is  based.  Also  the  con- 
vergence of  the  approximating  series  as  N -»  «>  appears  to  occur  in  an  irregu- 
lar or  oscillatory  fashion,  as  shown  by  the  three-mode  series  which  yields 
nearly  the  same  growth  curve  as  the  six-mode  series. 

Table  2.  Response  Functions  for  the  First  Six  Axial  Modes 


Mode 

IL 

3.596 

-1.221 

2L 

0.779 

-1.407 

3L 

0.459 

-.934 

4L 

0.357 

-.718 

5L 

0.306 

-.594 

6L 

0.273 

-.512 

A comparison  of  the  pressure  waveforms  obtained  with  the  four,  five, 
and  six  mode  expansions  is  shown  in  Figure  11.  The  principal  effect  of 
the  higher  harmonics  is  apparent  in  the  secondary  "wiggles"  which  are  depen- 
dent upon  the  number  of  modes  used  to  obtain  the  solution.  The  overall 
shape  of  the  pressure  waveforms  (i.e.,  the  steep  rise  to  maximum  and  more 
gradual  decline  to  minimum), the  frequency,  and  the  amplitude  are  less  strong- 
ly affected  by  the  number  of  modes. 

Figures  10  and  11  indicate  that  for  motors  without  particles,  more 
than  six  modes  (or  possibly  as  few  as  three  modes)  may  be  necessary  to  ade- 
quately represent  the  solution.  The  effect  of  the  number  of  modes  used 
to  describe  the  solution  for  cases  in  which  particles  are  present  will  be 
discussed  later  in  this  section. 

Effect  of  Initial  Disturbance.  This  study  was  designed  to  determine 
the  effect  of  the  amplitude  of  the  initial  disturbance  and  its  harmonic 
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content  upon  the  resulting  solutions.  Of  particular  interest  are  cases 
in  which  growth  to  limit-cycle  amplitude  occurs.  It  has  been  shown  in  Re- 
ferences (6),  (7),  and  (8)  for  liquid  rockets  that  the  final  limiting  ampli- 
tude is  independent  of  the  form  of  the  initial  disturbance.  On  the  other 
hand,  results  presented  by  Levine  and  Culick^  for  solid  rockets  indicate 
a significant  dependence  of  limit-cycle  amplitude  upon  initial  disturbance 
amplitude  for  amplitudes  in  the  range  5%  - 20%.  Thus  additional  data  is 
needed  to  clarify  this  issue. 

To  determine  the  effect  of  initial  disturbance  amplitude  upon  the 
resulting  limit-cycle  solution,  the  case  considered  previously . (no  particles, 

A = 6.00,  B = 0.590,  n = 0.583,  Q = 4.907)  was  run  for  initial  distur- 
bances of  3%,  7%,  and  20%.  Plots  of  head-end  pressure  amplitude  versus 
time  for  these  cases  are  shown  in  Figure  12  (using  a five  mode  series). 

The  3%  and  7%  solutions  appear  to  smoothly  approach  a limiting  amplitude 
of  12.5%  after  about  36  cycles.  The  20%  disturbance  also  approaches  the 
same  limiting  amplitude  but  in  an  oscillatory  fashion.  Although  the  paths 
are  different,  it  appears  that  if  the  calculations  are  continued  for  a 

sufficiently  long  time,  all  three  solutions  approach  the  same  limiting 

6 7 8 

amplitude.  This  result  is  consistent  with  previous  results  ’ ’ for  liquid 
rockets,  in  which  limit-cycle  amplitude  is  independent  of  the  initial  dis- 
turbance . 

For  a similar  case  for  which  2.5p,  particles  are  present,  3%  and 
10%  initial  amplitudes  yield  3.1%  and  9.0%  amplitude  oscillations  after 
12  cycles.  It  is  difficult  to  determine  whether  this  indicates  a true  depen- 
dence of  limit-cycle  amplitude  upon  initial  amplitude,  or  whether  the  approach 
to  limit-cycle  amplitude  is  just  extremely  slow.  Such  a slow  approach  to 
limiting  amplitude  is  expected  for  cases  in  which  the  linear  losses  due 
to  the  mean  flow,  particles,  and  the  nozzle  are  in  nearly  exact  balance 
with  the  gain  due  to  combustion  driving.  Such  a motor  is  said  to  be  operating 
near  the  point  of  linear  neutral  stability. 

The  effect  of  harmonic  content  of  the  initial  disturbance  was  studied 
by  introducing  initial  amplitudes  of  two  or  more  modes  simultaneously. 

Figure  13  shows  the  pressure  waveform  obtained  for  an  initial  disturbance 
composed  of  a combination  of  IL  and  2L  modes  for  the  case  of  no  particles. 

Due  to  nonlinear  coupling  of  modes,  the  amplitude  of  the  IL  mode  increases 
while  the  amplitude  of  the  2L  mode  decreases.  As  time  progresses  the  wave- 
form assumes  more  of  a IL  mode  character. 
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Number  of  Cycles 


Figure  12.  Effect  of  Initial  Disturbance  Amplitude  Upon  Approach 
to  Limiting  Amplitude 


Dimensionless  Time 


Effect  of  Propellant  Response  Function.  The  effect  of  propellant 
response  function  upon  the  stability  of  solid  rocket  motors  was  investigated 
using  both  approximate  and  "exact"  methods.  The  most  important  parameter 
influencing  the  amount  of  combustion  driving  was  found  to  be  the  real  part 
of  the  response  function  for  the  fundamental  mode,  This  depends  on 

the  values  of  A,  B,  and  n for  the  propellant,  and  the  frequency  parameter 
n.  Values  of  A,  B,  and  n were  selected  to  give  values  of  in  the  range 

0 to  4.  In  this  study  the  effect  of  particles  was  neglected. 

In  order  to  be  certain  that  the  propellant  response  function  is 

the  same  for  both  approximate  and  "exact"  solutions,  the  following  procedure 

* 

is  used.  First  a value  of  is  selected  for  input  into  the  "exact"  analy- 

sis; the  remaining  propellant  and  gas  properties  are  those  given  in  Table 
1.  The  corresponding  values  of  A,  B,  and  n are  then  calculated.  The  "exact" 
solutions  are  then  obtained  for  a given  initial  disturbance  and  the  resulting 
frequency  is  computed  from  the  pressure  waveforms  after  several  cycles. 

From  this  frequency  and  the  steady  state  burning  rate,  the  frequency  parameter 
n is  obtained  for  input  into  the  approximate  analysis.  The  approximate 
solutions  are  then  obtained  for  the  calculated  values  of  A,  B,  n,  and  for 
the  same  initial  disturbance  amplitude,  and  they  are  then  compared  with 
the  corresponding  "exact"  solutions. 

The  first  case  considered  was  a hypothetical  motor  in  which  the 
propellant  is  insensitive  to  the  pressure  oscillations;  that  is  R = 0 
for  all  modes.  Steady-state  combustion  was  included,  however,  giving  acous- 
tic losses  due  to  mean  flow,  flow  turning,  and  the  quasi-steady  nozzle. 

Exact  and  approximate  solutions  for  the  decay  of  the  pressure  amplitude 
with  time  are  shown  in  Figure  14  for  = 0.0847.  Excellent  agreement  be- 
tween the  approximate  and  "exact"  solutions  was  obtained  for  both  small 
(2.9%)  and  moderate  (10.4%)  amplitude  initial  di.sturbances . It  can  be  shown 
analytically  that  the  damping  rate  (.dimensionless)  approaches  a limiting 
value  at  small  amplitudes  given  by: 


O' 


(102) 


For  this  case  (t  = 1070  Hz)  the  "exact"  and  approximate  solutions  both  yield 
a damping  rate  of  about  -308  sec 
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Dimensionless  Time 

Figure  14.  Decay  of  Oscillations  Due  to  Mean  Flow,  Flow 
Turring,  and  Nozzle 
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In  the  second  case,  a value  of  = -131.8  cal/gm  was  assumed  yield- 
ing A = 6.00,  B = 0.643,  and  n = 0.616.  The  exact  calculations  gave 
= 4.91  for  which  “ 2.73.  Plots  of  pressure  amplitude  versus  time  for 
a 3.0%  initial  disturbance  is  shown  in  Figure  15  where  = 0.0780.  The 
initial  vertical  displacement  between  the  approximate  and  "exact"  solutions 
results  from  a mean  pressure  shift  obtained  with  the  "exact"  analysis  during 
the  first  cycle  of  oscillation  and  is  not  significant.  The  decay  rate  is 
much  less  than  in  the  previous  case  due  to  the  effect  of  combustion  driving; 
the  approximate  analysis  yields  a decay  rate  of  -14  sec  while  the  "exact" 
analysis  gives  -31  sec  Considering  that  the  net  damping  in  the  motor 

is  the  difference  between  two  relatively  large  quantities,  the  "exact"  and 
approximate  solutions  are  in  quite  good  agreement.  In  this  case  the  dis- 
crepancy in  the  damping  predicted  by  the  two  theories  is  17  sec  which 
is  only  about  6%  of  the  damping  in  the  absence  of  combustion  driving  (-288 
sec  ^ for  M = 0.078). 

e * 

Increasing  to  -134.7  cal/gm  gave  A = 6.00,  B = 0.607,  and  n = 
0.594.  For  a 3%  initial  disturbance  the  "exact"  analysis  gave  a growth 
to  a limiting  amplitude  of  about  2.8%  with  a corresponding  = 4.888  and 

R = 3.26.  For  the  approximate  analysis,  however,  the  oscillations  grew 

^ -1 
to  4.5%  after  12  cycles  and  were  still  growing  at  the  rate  of  36  sec 

For  this  case  solutions  for  a 10%  initial  disturbance  were  also  calculated; 
here  the  frequency  was  slightly  higher  giving  = 4.932  and  = 3.23. 

Both  solutions  were  still  decaying  after  12  cycles,  the  "exact"  solution 
at  7.5%  amplitude  with  » = -42  sec  and  the  approximate  solution  at  9.5% 
amplitude  with  ot  = -31  sec  A limit  cycle  of  about  7%  for  the  approxi- 
mate analysis  is  consistent  with  these  results. 

★ 

Finally,  a value  of  = -136.1  cal/gm  was  chosen  for  which  A = 

6.00,  B = 0.590,  and  n=  0.583.  Amplitude-time  plots  are  given  for  3%  and 
7%  initial  amplitudes  in  Figure  16,  for  which  = 3.6.  This  figure  shows 
that  again  the  approximate  analysis  predicts  a somewhat  greater  initial 
growth  rate  and  larger  limit-cycle  amplitude  than  the  "exact"  analysis. 

In  this  case  the  approximate  analysis  yields  a growth  rate  of  73  sec  ^ for 
the  3%  solutions  and  a final  amplitude  of  11.3%,  while  the  corresponding 
values  obtained  with  the  "exact"  analysis  are  38  sec  ^ and  8.0%  respectively. 
The  frequencies  obtained  are  in  close  agreement  at  about  1050  Hz.  Head- 
end  and  mid-chamber  waveforms  are  shown  in  Figure  17  for  limit-cycle  con- 


70 


1.0 

0.8  - 

0.6  - 

0.4  - 

0.2  - 

0.1. 

i-  0.08. 


0) 

T3 

P 

4J 


(1) 

>-l 

P 

W 

CO 

0) 

Ui 

•o 

c 

UJ 

I 

•a 

<D 

X 

jti 

(U 


1 

cd 

CJ 

cu 


Galerkin  Method  (6  Modes) 


Exact  Analysis 


3%  Initial  Amplitude 


0.06U. 


0.04L 


0.02L 


0.01. 

0.008. 


0.006L 


0.004L 


0.002h 


No  Particles 

A = 6.00 

B = 0,643 

n = 0.616 

0 = 4.91 

\ = 1.23 

M = 0.078 
e 


O.OOlL 


0 


Figure  15. 


8 12  16 
Dimensionless  Time 


20 


24 


Decay  of  Oscillations  in  a Stable  Motor  Without 
Particles 


71 


Pressure  Perturbation 


0.1 


Figure  17.  Pressure  Waveforms  at  Limiting  Amplitude  for  a 
Motor  Without  Particles 


ditions  (after  11  cycles  of  oscillation).  The  approximate  waveforms  show 
considerably  more  harmonic  content  than  the  "exact"  waveforms,  but  the  quali- 
tative features  of  the  waveforms  are  in  good  agreement  (i.e.,  waveform 
steepening  and  phase  relation  between  even  harmonics  and  the  fundamental). 

The  Effect  of  Oscillation  Frequency.  The  dependence  of  unstable 
solid  rocket  limit-cycle  behavior  upon  the  oscillation  frequency  was  studied 
by  "arying  the  frequency  parameter  Q for  a given  propellant.  Figures  18 
and  19  describe  the  head-end  pressure  oscillations  for  two  rockets  having 
identical  propellants  (A  = 7.54,  B = 0.880,  n = 0.85)  but  different  nondimen- 
sional  frequencies  Cl.  The  magnitudes  of  the  combustion  response  "driving" 
for  the  five  modes  present  in  these  two  cases  are  also  included  in  these 
figures  which  clearly  show  that  for  the  case  = 1.64  (Figure  18)  the 
magnitudes  of  the  combustion  responses  of  the  higher  harmonics  are  considerably 
larger  than  in  the  case  when  0^^  = 10.0  (Figure  19).  This  difference  in 
driving  is  directly  responsible  for  the  observed  difference  in  the  head- 
end  pressure  wave  forms.  Furthermore,  these  figures  suggest  that  more  than 
five  modes  are  probably  needed  to  adequately  describe  the  solution  in  the 
= 1.64  case. 

Effect  of  Particles.  In  the  last  phase  of  the  parametric  study 
the  effect  of  particle  size  and  concentration  upon  the  nonlinear  stability 
characteristics  of  solid  rocket  motors  was  investigated.  Three  basic  cases 
were  considered:  (1)  attenuation  of  waves  in  a particle-gas  mixture  in 

a closed-ended  chamber  (i.e.,  the  case  of  particles  in  a box  treated  by 
Culick^^  and  Levine  and  Culick^),  (2)  attenuation  of  waves  in  a motor  with- 
out combustion  driving  ( R.  = 0),  and  (3)  limit-cycle  solutions  in  an  un- 
stable motor.  In  all  cases  approximate  solutions  obtained  with  the  Galerkin 
method  (six  modes)  are  compared  with  "exact"  solutions  obtained  with  the 
Kooker-Zinn  model.  For  the  case  of  particles  in  a box,  solutions  obtained 
with  Culick's  MOA  and  Levine's  "exact"  analysis  are  also  included. 

For  the  case  of  particles  in  a box,  the  input  parameters  are  the 
same  as  those  given  by  Culick.^^  These  are  needed  in  order  to  calculate 
the  particle  drag  constant  K for  use  in  the  "exact"  and  approximate  analyses. 
Since  Culick  specifies  the  equilibrium  frequency  for  the  particle-gas  mix- 
ture rather  than  the  chamber  length  and  speed  of  sound,  the  value  of  K for 
a given  particle  size  was  calculated  from  the  following  expression: 
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Head-End  Pressure  Perturbation 


(103) 


K = 


£* 

g 


where  /2L  is  the  frequency  for  the  pure  gas.  The  gas  viscosity 

was  obtained  from  the  formula  given  in  Reference  10: 


kg/m-sec 


(104) 


for  T = 3416  K used  by  Culick.  In  the  approximate  analysis  the  ratio  of 
specific  heats  Y for  the  mixture  replaces  Y for  the  gas.  The  value  of 
Y depends  on  the  specific  beats  of  the  gas  and  particle  material  and  the 
concentration  of  the  particles  and  is  given  by: 


Y 


(105) 


where  C and  are  the  specific  heats  of  the  particles  and  gas  respectively. 
Also  in  specifying  the  initial  disturbance,  the  equilibrium  frequency  is 
needed;  it  is  related  to  the  pure  gas  frequency  by: 


f = 
e 


(106) 


To  evaluate  the  accuracy  of  the  approximate  analysis  in  describing 
the  effects  of  linear  particle  damping,  "exact''  and  approximate  solutions 
were  obtained  for  the  case  of  a particle-gas  mixture  in  a closed-ended  box. 
Assuming  an  initial  disturbance  amplitude  of  3Z,  values  of  damping  and  frequency 
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The  results 


were  computed  for  several  values  of  the  particle  diameter  a . 
were  calculated  for  a pure-gas  frequency  of  f = 1071  Hz  and  a particle 
concentration  of  = p^/p  = 0.20.  The  values  of  K corresponding  to  the 
particle  sizes  considered  in  this  study  are  given  in  Table  3. 

Table  3.  Drag  Constants  for  Approximate  and  "Exact"  Models 


Particle  Size,  a 

K 

k' 

(p.) 

(Approximate) 

(Exact) 

o 

• 

CM 

46.72 

160.2 

2.5 

29.90 

102.5 

3.0 

20.76 

71.18 

4.0 

11.68 

40.04 

6.0 

5.191 

17.80 

8.0 

2.920 

10.01 

9.0 

2.307 

7.909 

10.0 

1.869 

6.406 

15.0 

.8306 

2.847 

20.0 

.4672 

1.602 

25.0 

.2990 

1.025 

30.0 

.2076 

.7118 

40.0 

.1168 

.4004 

The  value  of  used  in  the  "exact"  analysis  differs  from  the  value  K 
used  in  the  approximate  analysis  due  to  the  difference  in  the  reference 
state  used  in  nondimensionalizing  the  "exact"  equations.  In  the  Kooker- 
Zinn  analysis  standard  atmospheric  conditions  are  used  as  the  reference 
state,  while  the  chamber  stagnation  conditions  are  used  in  the  approximate 
analysis.  Thus  the  two  particle  drag  constants  are  related  by 

k'  = K (c*/c*  ,)  (107) 

o ref  ' ' 

* ★ 

where  c and  c ^ are  the  speeds  of  sound  at  chamber  stagnation  condi- 
o ref  ® 

tions  and  standard  atmospheric  conditions  respectively. 

Curves  of  decay  rate  and  frequency  as  a function  of  particle  diameter 
are  presented  in  Figure  20  for  particles  in  a box  (l.e.,  = 0),  The 

damping  and  frequency  shift  arise  solely  from  the  particle-gas  interaction. 


( 

1 
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Decay  Rate,  -a(sec 


Frequency,  f (Hz) 


since  mean  flow,  combustion,  and  nozzle  effects  are  absent.  The  damping 
is  low  for  particles  smaller  chan  2p,  and  for  particles  larger  than  30 
Pi  , and  the  maximum  damping  occurs  at  about  8.5p,  for  the  frequency  con- 
sidered. The  frequency  varies  from  the  equilibrium  value  of  985  Hz  for 
small  particles  ( a <4pi  ) to  the  pure-gas  value  of  1080  Hz  ( > 25p>  ). 

This  figure  shows  very  good  agreement  between  the  decay  rates  and  frequencies 
predicted  by  the  approximate  and  "exact"  models  over  the  size  range  from  2 p, 
to  40  p . This  indicates  that  the  approximate  analysis  correctly  models 
the  particle-gas  interaction  for  small  amplitude  disturbances  in  the  absence 
of  mean  flow  and  combustion  effects. 

To  further  investigate  the  particle-gas  interaction,  several  of 

the  cases  considered  by  Culick^^  were  calculated  using  the  Galerkin  Method 

and  the  Kooker-Zinn  "exact"  analysis.  These  solutions  were  then  compared 

with  Culick's  MOA  solutions  and  Levine's  "exact"  solutions  for  the  same 

cases.  The  case  of  2.5  p particles  with  C = 0.36  and  f = 800  Hz  was 

me 

first  considered  (K  = 32.9).  Plots  of  damping  versus  number  of  cycles 
for  3%  initial  disturbances  are  presented  in  Figure  21,  while  similar  plots 
for  15%  initial  disturbances  are  shown  in  Figure  22.  Figure  21  shows  very 
good  agreement  between  the  Galerkin  solutions  and  the  Kooker-Zinn  solutions, 
while  the  agreement  beween  the  Galerkin  method  and  Culick's  and  Levine's 
solutions  is  not  as  good.  For  15%  initial  disturbances.  Figure  22  shows 
that  the  Galerkin  solutions  approach  the  "exact"  solutions  asymptotically, 
although  significant  quantitative  differences  in  the  computed  decay  rates 
occur  during  the  initial  cycles  of  oscillation.  This  figure  also  shows 
the  effect  of  number  of  modes  upon  the  computed  decay  rates;  in  this  case 
the  solution  is  adequately  described  by  a five-mode  series.  The  Galerkin 
solutions  for  the  15%  disturbance  are  also  compared  with  Culick's  (MOA) 
and  Levine's  ("exact")  solutions  in  Figure  22.  Again  the  agreement  with 
Culick's  and  Levine's  solutions  is  not  as  good  as  with  Kooker ' s solutions. 

Some  of  the  discrepancies  between  the  Galerkin  solutions  and  Culick's 
MOA  solutions  (as  well  as  Levine's  "exact"  solutions)  may  be  due  to  differences 
in  the  physical  models  involved.  For  example,  the  Galerkin  method  results 
shown  were  obtained  using  a linear  drag  law  to  describe  the  particle  damping, 
while  the  results  of  Culick  and  Levine  were  obtained  using  a nonlinear  drag 
law.  Furthermore,  the  analyses  of  Culick  and  Levine  included  a thermal 
equation  for  the  particles,  while  this  equation  was  ignored  in  both  the 
Galerkin  and  the  Kooker-Zinn  analyses. 
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Two  additional  cases  were  considered  for  particles  in  a box:  (1) 

a = 2.5p,  , fg  = 1500  Hz  (K  = 17.55)  and  (2)  a = 10  p,  , = 1500  Hz  (K 

= 1.097).  For  both  cases  C = 0.36  and  a 15%  initial  disturbance  was  consid- 

m 

ered.  Approximate  and  "exact"  solutions  for  these  cases  are  shown  in  Figure 

23.  For  the  2.5 p particles  the  agreement  between  approximate  and  "exact" 

solutions  at  1500  Hz  is  better  than  at  800  Hz  (Figure  22),  and  a four-mode 

series  is  adequate  to  describe  the  solution.  For  the  10 p case  both  theories 

predict  large  variations  in  decay  rate  from  cycle  to  cycle,  and  significant 

differences  between  four,  five,  and  six  mode  solutions  are  apparent.  Although 

for  any  given  cycle  the  decay  rates  predicted  by  the  approximate  and  "exact" 

analyses  can  differ  significantly,  the  mean  decay  rates  during  the  first 

14  cycles  are  in  fairly  good  agreement. 

The  second  case  considered  in  the  study  of  particle  damping  effects 

was  a motor  for  which  R = 0,  that  is  the  combustion  process  does  not  respond 

to  the  pressure  oscillations.  Thus  the  effects  of  mean  flow  and  nozzle 

damping  are  considered  along  with  the  particle  damping  effects.  Calculations 

of  decay  rate  and  frequency  were  made  for  a motor  with  = 0.0780  and  f^ 

= 1071  Hz  for  several  values  of  a and  C . 

m 

The  effect  of  particle  size  on  the  "exact"  and  approximate  solutions 

for  C = 0.20  is  shown  in  Figure  24  for  3%  initial  disturbances.  These 
m 

curves  are  similar  to  those  shown  in  Figure  20  for  particles  in  a box,  except 
that  the  damping  curves  are  shifted  upwards  due  to  the  large  damping  due 
to  mean  flow  and  nozzle  effects  i a = -289  sec  ^ in  the  absence  of  particles). 
The  decay  rate  curves  show  that  the  approximate  and  "exact"  analyses  are 
in  excellent  agreement  for  particle  diameters  larger  than  about  15p,  , and 
both  analyses  predict  practically  the  same  optimum  particle  size  for  maximum 
damping.  For  particles  smaller  than  15|i  the  approximate  analysis  predicts 
a smaller  linear  decay  rate  than  the  "exact"  analysis  with  a maximum  error 
of  about  13%.  This  discrepancy  is  probably  due  to  certain  terms  proportional 
to  the  particle  drag  constant  K and  the  mean  flow  Mach  number  (see  Appendix 
A)  that  were  neglected  as  higher  order  in  the  approximate  analysis.  Since 
K is  inversely  proportional  to  the  frequency  of  oscillation  and  the  square 
of  the  particle  diameter,  better  agreement  is  obtained  in  the  high-frequency 
and  large-particle  size  range.  The  frequencies  predicted  by  the  two  models 
are  in  best  agreement  for  small  ( a < 6 U ) and  large  particles  ( a > 20p,  ). 

It  should  be  noted  that  the  good  agreement  between  "exact"  and  appro- 
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Figure  24.  Effect  of  Particle  Size  on  Decay  Rate  and  Frequency 
for  Motor  Without  Combustion  Driving 
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ximaCe  analyses  for  large  particles  has  been  obtained  under  the  assumption 
that  the  particle-gas  interaction  is  described  by  the  Stokes  Drag  Law  (Equa- 
tions (7)  and  (8)).  For  large  particle  sizes  and  high  frequencies  (i.e., 
large  Reynolds  numbers),  however,  the  Stokes  Drag  Law  is  no  longer  valid 
and  the  nonlinear  drag  law  given  by  Equation  (79)  should  be  used.  The  case 
of  nonlinear  particle  damping  is  considered  in  a later  subsection. 

The  dependence  of  damping  and  frequency  upon  the  particle  concentra- 
tion C is  shown  in  Figure  25  for  10  p,  particles.  Both  approximate  and 

ID 

"exact"  analyses  predict  a nearly  linear  dependence  of  frequency  and  damping 
upon  particle  concentration  for  small  amplitude  initial  disturbances  (3Z). 

Here  the  decay  rate  increases  from  -289  sec  ^ for  no  particles  (C  = 0) 

-1  ™ 
to  about  -800  sec  for  = 0.35,  while  the  frequency  decreases  from  the 

pure-gas  value  of  1071  Hz  to  about  1010  Hz  at  = 0.35.  For  this  particle 
size,  the  approximate  analysis  yields  smaller  decay  rates  than  the  "exact" 
analysis,  with  a discrepancy  of  about  5%  at  = 0.35.  The  predicted  fre- 
quencies are  also  lower  for  the  approximate  analysis,  but  the  discrepancies 
are  all  less  than  1%. 

The  last  case  considered  in  the  particle  damping  study  was  a motor 

with  linear  combustion  driving.  The  same  basic  motor  parameters  were  used 

as  for  the  cases  considered  previously  when  particles  were  absent.  As  in 

Levine's  cases^,  a particle  diameter  of  2.5(1  was  selected  at  a concentration 

of  0.36.  In  order  to  obtain  a limiting  amplitude,  a larger  propellant  response 

than  that  of  Figure  16  was  needed  to  overcome  the  increased  damping  due 

★ 

to  the  particles.  Thus  a value  of  = -137.0  cal/gm  was  chosen  which  gave 
A = 6.00,  B = 0.580,  n = 0.575.  The  corresponding  response  curve  has  a 
peak  of  = 4.2  at  Q =4.5. 

Plots  of  head-end  pressure  amplitude  versus  time  for  this  case  are 
presented  in  Figure  26  for  3%  and  10%  initial  disturbances.  The  agreement 
between  the  approximate  and  "exact"  analyses  is  seen  to  be  fairly  good. 

Both  analyses  yield  a limit-cycle  of  about  9%  amplitude  for  the  10%  initial 
disturbance  and  a relatively  slow  growth  for  the  3%  initial  disturbance. 

For  the  low  amplitude  disturbance,  the  approximate  analysis  gives  the  smaller 
growth  rate,  which  is  in  contrast  to  the  case  when  particles  are  absent 
(Figure  16)  where  the  approximate  analysis  predicted  a larger  growth  rate 
than  the  "exact"  analysis.  This  result  was  unexpected,  since  the  approxi- 
mate analysis  underestimates  the  damping  for  the  case  of  2 . 5 p.  particles 
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without  combustion  driving  (Figure  24).  The  corresponding  head-end  and 
mid-chamber  pressure  waveforms  for  limit-cycle  conditions  are  shown  in 
Figure  27.  Both  analyses  yield  nearly  sinusoidal  waveforms;  the  only  evi- 
dence of  higher  harmonics  is  the  steepening  of  the  rising  branch  of  the 
waveform  and  the  double  frequency  oscillation  at  the  center  of  the  chamber. 

This  is  due  to  the  presence  of  particles  which  suppress  the  higher  harmonics. 

The  importance  of  using  the  proper  frequency  parameter  in  the  appro- 
ximate analysis  when  making  comparisons  with  the  "exact"  analysis  is  readily 
illustrated  by  this  case.  From  the  "exact"  analysis,  the  frequency  parameter 
for  the  10%  initial  amplitude  is  0 - 4.244  for  which  ^ “ 4.17.  If, 

however,  the  pure-gas  frequency  is  assumed,  ” 4.94  for  which  = 3.78, 
and  the  approximate  analysis  would  predict  a decaying  oscillation.  Thus 
the  approximate  analysis  would  appear  to  give  the  wrong  results,  because 
the  response  function  was  not  the  same  in  both  cases.  For  this  reason  the 
"exact”  program  must  always  be  run  first  to  obtain  the  proper  Q for  input 
into  the  approximate  program.  Of  course,  if  the  calculations  are  being 
made  with  the  approximate  analysis  only,  the  values  of  A,  B,  n,  and  0 
must  be  estimated  or  the  values  of  R for  the  various  modes  may  be  obtained 
from  experimental  data. 

Some  parametric  variations  on  the  case  considered  above  will  now 
be  considered.  The  effect  of  variations  of  particle  size  upon  the  limit- 
cycle  amplitude  and  waveform  are  shown  in  Figure  28  for  = 0.10.  Here 
head-end  pressure  waveforms  are  shown  for  a - 2.5,  3.5,  and  5 ; the  corres- 
ponding amplitudes  are  19.3%,  15.3%,  and  8.2%.  This  decrease  in  limiting 
amplitude  with  increasing  particle  size  is  consistent  with  the  increasing 
particle  damping  with  increasing  tj  as  shown  in  Figure  20.  The  distortion 
of  the  waveform  also  decreases  with  increasing  a due  to  the  decreasing 
amplitude  and  the  increased  attenuation  of  the  higher  harmonics  by  the 
particles.  For  particles  with  diameters  somewhat  greater  than  5^  (for  the 
given  propellant  response  and  particle  concentration),  limit-cycles  are 
not  obtained  and  the  oscillations  decay.  For  continued  increases  in  parti- 
cle size  the  decay  rate  reaches  a maximum,  then  declines,  and  limit-cycles 
are  again  obtained  for  large  particles  (a  greater  than  roughly  15|j,  ).  For 
large  particles,  the  opposite  trend  is  obtained,  that  is,  limit-cycle  ampli- 
tude increases  with  increasing  particle  size. 

The  influence  of  particle  concentration  C upon  the  limit-cycle 


L 
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amplitude  and  pressure  waveform  is  shown  in  Figures  29  and  30  for  2.5p,  parti- 
cles. Figure  29  gives  limit-cycle  pressure  waveforms  for  C = 0.10,  0.15, 

m 

and  0.20  for  which  the  corresponding  amplitudes  are  19.3%,  16.5%,  and  14.3% 
respectively.  The  increased  attenuation  of  the  higher  harmonics  with  increas- 
ing particle  concentration  is  evident  from  these  waveforms.  A plot  of 
limiting  amplitude  as  a function  of  particle  concentration  is  presented 

in  Figure  30  for  0.05  ^ C S 0.03.  For  C >0.1  the  decline  in  amplitude 
° m m 

with  increasing  C is  seen  to  be  nearly  linear,  while  a sharp  rise  in 
m 

amplitude  occurs  for  C <0.1. 

^ m 

In  studying  the  effe<'ts  of  a and  C upon  limiting  amplitude  and 

m 

waveform,  the  frequency  parameter  was  held  constant  at  Q = 4 . 2 in  order 

to  isolate  the  effects  of  particle  damping.  In  an  actual  motor  the  situation 

is  more  complicated.  The  oscillation  frequency  is  dependent  on  particle 

size  and  concentration  (Figures  24  and  25),  thus  changing  a or  C also 

m 

changes  the  propellant  response  function  A through  the  effect  of  frequency 
on  Q . Depending  on  the  values  of  A,  B,  and  n for  the  given  propellant, 
allowing  0 to  vary  will  yield  higher  or  lower  limit-cycle  amplitudes  than 
would  be  obtained  for  * constant.  It  should  be  noted  from  Figure  24 

that  for  small  and  large  particles  the  equilibrium  and  pure-gas  frequencies, 
respectively,  are  good  approximations  for  use  in  calculating  0 . 

Sumnary  of  Mode-Coupling  Results.  The  principal  results  of  this 
study  to  assess  the  importance  of  gasdynamic  mode-coupling  upon  axial  mode 
instabilities  in  solid  rocket  motors  will  now  be  summarized.  It  has  been 
shown  that  the  solutions  obtained  with  the  approximate  and  "exact"  analyses 
are  in  fairly  good  agreement  regarding  growth/decay  rates,  frequencies, 
limiting  amplitudes  and  pressure  waveforms.  These  results  indicate  that 
the  approximate  analysis  (Galerkin  Method)  accounts  for  mean  flow,  flow- 
turning, and  nozzle  damping  effects,  as  well  as  the  effects  of  particle 
damping,  combustion  driving,  and  nonlinear  gasdynamic  mode-coupling.  These 
studies  have  shown  that  gasdynami:  mode-coupling  is  probably  the  most  impor- 
tant nonlinear  process  which  influences  unstable  motor  behavior.  Gasdynamic 
mode-coupling  appears  to  limit  the  oscillation  amplitude  in  unstable  motors 
by  the  transfer  of  energy  from  the  unstable  fundamental  mode  to  the  stable 
higher  frequency  modes.  On  the  other  hand,  the  lack  of  pulsed  instabilities 
in  the  results  presented  in  this  section  indicates  that  gasdynamic  nonlinear- 
ities (at  least  to  second  order)  can  not  account  for  pulsed  or  triggered 
instabilities.  These  observations  on  the  role  of  gasdynamic  mode-coupling 
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in  solid  rocket  motors  are  in  agreement  with  the  results  of  previous  studies 
for  liquid  rockets^’^. 

Concerning  the  effect  of  number  of  modes  upon  the  approximate  solu- 
tions, the  following  results  are  notable.  The  number  of  modes  necessary 
to  adequately  represent  the  solution  depends  upon  the  relative  stability 
of  the  higher  frequency  modes  with  respect  to  the  fundamental  mode.  For 
cases  in  which  the  higher  modes  are  heavily  damped,  as  in  the  case  when 
small  particles  are  present,  a three  or  four  mode  series  expansion  appears 
to  be  aot.quate  (Figures  22  and  23).  For  motors  with  large  particles  or 
for  motors  without  particles,  the  higher  modes  are  less  strongly  damped 
and  six  or  more  modes  may  be  required  (Figures  10  and  23).  For  cases  in 
which  the  propellant  response  for  the  higher  modes  are  comparable  to  or 
greater  than  the  propellant  response  for  the  fundamental  mode  (Figure  18), 
the  convergence  of  the  series  expansion  may  be  so  slow  that  an  excessive 
number  of  modes  are  required.  In  the  latter  case,  the  approximate  technique 
may  require  more  computation  time  than  the  "exact"  analysis. 

Other  results  of  the  parametric  studies  are  of  interest.  For  the 
cases  considered,  limiting  amplitude  appeared  to  be  independent  of  initial 
disturbance  amplitude  and  harmonic  content.  The  linear  combustion  response 
factor  is  a major  factor  influencing  growth/decay  rates  and  limiting 

amplitudes.  For  cases  without  particles,  oscillations  decay  for  < 3.3 

while  for  >3.3  limiting  amplitude  increases  with  increasing  R^  . 

For  cases  with  particles,  larger  values  of  R^  are  needed  for  instability. 

The  oscillation  frequency,  determined  by  the  chamber  length  and  sound  speed, 
affects  the  stability  characteristics  of  the  motor  primarily  through  its 
influence  on  combustion  driving  and  particle  damping.  Finally,  particle 
size  and  concentration  are  both  important  parameters  influencing  growth  and 
decay  rates,  limiting  amplitudes  and  waveforms;  there  is  an  optimum  value 
of  the  particle  drag  constant  K (determined  by  particle  size  and  frequency) 
for  maximum  damping. 

4 . 3 Effect  of  Nonlinear  Combustion  Driving 

In  the  studies  presented  so  far,  it  was  assumed  that  the  only  non- 
linearity present  in  the  solid  rocket  system  was  gasdynamic  mode-coupling. 

In  the  next  phase  of  this  investigation,  the  effects  of  additional  nonlinear- 
ities were  included  in  the  analysis.  The  effects  of  nonlinearities  in 
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the  pressure-coupled  transient  combustion  process  upon  the  stability  charac- 
teristics of  solid  rocket  motors  are  presented  and  discussed  in  this  section. 
Velocity-coupled  nonlinear  combustion  processes,  which  are  also  known  to 
be  important,  were  not  considered  in  this  study  due  to  the  time  and  economic 
limitations  of  this  project. 

The  results  presented  in  this  section  were  obtained  by  the  "exact" 
Kooker-Zinn  model  using  the  nonlinear  combustion  option  and  by  the  approxi- 
mate analysis  (Galerkin  method)  using  the  heuristic  nonlinear  combustion 

3 5 

model  described  in  Section  2.7.  Kooker  and  Zinn  and  Levine  and  Culick 
have  boi,h  considered  the  effects  of  combustion  nonlinearities  using  "exact" 
numerical  methods,  and  both  have  shown  that  the  effect  of  a nonlinear  com- 
bustion response  is  significant  for  amplitudes  greater  than  about  7 - 10%. 

On  the  other  hand,  this  study  is  the  first  attempt  at  incorporating  burning 
rate  nonlinearitieo  into  an  approximate  instability  model.  One  of  the  major 
objectives  of  this  study,  therefore,  is  to  assess  the  feasibility  of  modeling 
combustion  nonlinearities  using  the  approximate  technique. 

The  approximate  model  was  first  used  to  study  the  basic  case  considered 
previously  in  the  mode-coupling  investigations.  The  linear  combustion  parameters 
for  this  case  are  A = 6.00,  B = 0.59,  n = 0.583;  there  are  no  particles 
present  in  the  flow;  and  the  mean  flow  Mach  number  at  the  nozzle  entrance 
is  0.0780.  The  acoustic  frequency  and  the  steady-state  burn  rate  yield 
Q = 4.907.  The  nonlinear  combustion  parameter  b^  appearing  in  the  heuris- 

tic nonlinear  combustion  moael  (Equation  (77))  is  not  known  a priori.  In 
general,  b^  is  a complex  number  which  is  expected  to  be  different  for  each 
acoustic  mode.  In  order  to  reduce  the  number  of  parameters  to  be  considered 
in  this  study,  it  was  assumed  that  the  are  all  real  numbers  and  that 

they  are  the  same  for  all  modes.  Thus  the  effect  of  combustion  nonlinearities 
was  studied  by  varying  the  single  parameter  b in  the  approximate  analysis. 

Approximate  solutions  for  the  growth  of  a 7%  initial  disturbance 
to  limiting  amplitude  are  shown  in  Figure  31  for  both  linear  and  nonlinear 
combustion  models.  For  the  nonlinear  combustion  response,  both  positive 
(b  = 0.5)  and  negative  (b  = -0.5)  values  of  b were  considered.  This  figure 
shows  a significant  effect  of  combustion  n'^nlinearities  upon  the  initial 
growth  rate  and  limit-cycle  amplitude  of  the  pressure  oscillations  in  the 
motor.  For  the  positive  value  of  b,  both  the  initial  growth  rate  and  limit- 
cycle  amplitude  are  greater  than  the  corresponding  values  obtained  with 
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' the  linear  combustion  response  (b  = 0),  while  the  negative  value  of  b yields 

i reduced  values  of  the  initial  growth  rate  and  limiting  amplitude.  The 

' limiting  amplitudes  obtained  were  as  follows:  11.3%  for  linear  combustion, 

10.1%  for  b = -0.5,  and  12.9%  for  b = 0.5.  This  behavior  is  expected  from 
I the  form  of  the  assumed  nonlinear  combustion  model  (Equation  (77)),  since  for 

b > 0 the  effective  combustion  response  increases  with  increasing  ampli- 

1 tude  and  for  b < 0 the  combustion  response  decreases  with  increasing  ampli- 

j tude. 

Figure  32  shows  the  dependence  of  limit-cycle  pressure  amplitude 
upon  the  nonlinear  combustion  parameter  b for  the  case  considered  above. 

For  -1.0  s b i 1.0  it  is  seen  that  limiting  amplitude  increases  nonlinearly 
with  increasing  b and  that  the  strongest  effect  of  nonlinear  combustion 
occurs  for  positive  values  of  b.  Two  factors  contribute  to  this  nonlinear 
dependence  of  limiting  amplitude  upon  b:  (1)  limiting  amplitude  is  a non- 
linear function  of  and  (2)  is  proportional  to  the  pressure  ampli- 

tude as  well  as  the  value  of  b. 

In  order  to  determine  realistic  values  of  b to  use  in  the  approxi- 
mate analysis,  comparisons  were  made  with  "exact"  solutions  obtained  with 

the  nonlinear  burning  rate.  "Exact"  solutions  obtained  with  linear  and 

★ 

nonlinear  burning  rates  are  compared  in  Figure  33  for  the  same  case  (Q^ 

= -136.1  cal/gm)  considered  above.  Here  it  is  seen  that  nonlinear  combus- 
tion effects  yield  a larger  limiting  amplitude  (8.7%  nonlinear  versus  7.9% 
linear).  It  should  be  noted  that  for  this  case  the  approximate  analysis 
yields  a significantly  larger  limiting  amplitude  (11.3%)  than  the  "exact" 
analysis  for  linear  combustion,  therefore  the  appropriate  value  of  b cannot 
i be  obtained  by  direct  comparison  of  approximate  and  "exact"  solutions. 

Instead,  a value  of  b is  chosen  which  gives  the  same  relative  increase  in 
limiting  amplitude  due  to  combustion  nonlinearities  as  the  "exact"  analysis. 

; For  this  case,  nonlinear  combustion  increases  the  limiting  amplitude  by 

t 

I a factor  of  1.10  according  to  the  "exact"  analysis.  Thus  the  corresponding 

‘ approximate  solution  should  yield  a limit-cycle  amplitude  of  12.4%  corres- 

ponding to  b = 0.32  (Figure  32). 

The  effects  of  a nonlinear  burning  rate  upon  the  pressure  and  burn- 
ing rate  waveforms  are  shown  in  Figure  34  for  the  "exact"  analysis,  while 
the  corresponding  approximate  solutions  for  the  pressure  waveforms  (b  = 

0.32)  are  given  in  Figure  35.  Figure  34  shows  that  a nonlinear  combustion 
response  modifies  significantly  the  burning  rate  waveforms;  the  positive 

98 


part  of  the  waveform  becomes  steeper  with  higher  peak  burning  rates  and 
the  negative  part  becomes  flatter.  On  the  other  hand,  both  approximate 
and  "exact"  models  predict  that,  for  the  moderate  amplitudes  obtained  in 
this  case,  combustion  nonlinearities  have  very  little  effect  on  the  shape 
of  the  pressure  wave  (Figures  34  and  35).  The  only  noticeable  effect  of 
nonlinear  burning  rate  upon  the  pressure  waveform  is  a small  mean  pressure 
shift  which  is  obtained  with  the  "exact"  analysis  but  not  with  the  approxi- 
mate model. 

The  steepening  of  the  burning  rate  waveforms  and  the  insensitivity 
of  the  pressure  wave  shape  to  combustion  nonlinearities  are  in  agreement 
with  previously  reported  "exact"  solutions  obtained  by  Levine  and  Culick.^ 
However,  Levine  and  Culick  considered  a different  set  of  combustion  parameters 
(A  = 6.0,  B = 0.53,  n = 0.3)  and  included  2 p.  particles  in  the  flow,  which 
gave  a reduction  in  limiting  amplitude  from  31%  to  15%  when  burning  rate 
nonlinearities  were  included  in  the  analysis.  This  trend  is  opposite  from 
that  obtained  above,  and  it  indicates  that  negative  values  of  b may  be 
appropriate  in  some  cases. 

To  further  investigate  the  effect  of  combustion  nonlinearities  on 
the  stability  of  motors,  a second  case  was  considered  which  was  obtained 
by  adding  2.5  p diameter  particles  (C^  = 0.36)  to  the  case  considered  above. 
The  linear  combustion  parameters  (A  = 6.00,  B “ 0.59,  n = 0.583)  and  the 
steady-state  Mach  number  at  the  nozzle  entrance  (M^  ” 0.0780)  remain  the 
same  as  before.  The  oscillation  frequency,  however,  is  lowered  from  1071 
Hz  to  918  Hz  due  to  the  effect  of  the  particles.  This  decrease  in  frequency 
results  in  an  increase  in  the  linear  combustion  response  factor  for  the 
fundamental  mode;  that  is,  increases  from  3.60  to  3.81.  The  increased 

damping  due  to  the  presence  of  the  particles  offsets  the  increased  combustion 
driving  to  yield  a system  which  is  stable  to  moderate  amplitude  disturbances. 
Therefore,  in  the  results  to  follow  the  effects  of  combustion  nonlinearities 
upon  the  decay  rate,  rather  than  limiting  amplitude,  are  presented. 

The  effect  of  the  burning  rate  nonlinearities  upon  the  decay  rates 
for  3%  and  10%  initial  disturbances  is  shown  in  Figure  36  as  calculated 
by  the  "exact"  analysis.  Here  the  nonlinear  combustion  response  yields 
larger  decay  rates,  which  reflects  a decrease  in  the  amount  of  combustion 
driving  with  the  addition  of  nonlinear  combustion  effects.  For  the  3%  dis- 
turbance, combustion  nonlinearities  increase  the  damping  after  10  cycles 
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from  -2.8  sec  ^ to  -11.1  sec  ^ (-8.3  sec  ^ increase),  while  an  increase 
in  damping  from  -16.7  sec  ^ to  -28.4  sec  ^ (-11.7  sec  ^ increase)  occurs 
for  the  10!{  disturbance.  This  trend  of  decreased  combustion  driving  with 
the  addition  of  combustion  nonlinearities  is  opposite  to  the  trend  obtained 
earlier  for  the  motor  without  particles  and  in  agreement  with  the  results 
of  Levine  and  Culick^  for  a motor  with  2 p,  - diameter  particles.  It  is 
believed  that  this  change  in  the  behavior  of  the  nonlinear  combustion  re- 
sponse when  particles  are  included  is  simply  a result  of  the  decreased  fre- 
quency of  oscillation,  which  apparently  affects  the  nonlinear  character- 
istics of  the  combustion  response  as  well  as  the  linear  burning  characteris- 
tics. These  results  also  show  that  combustion  nonlinearities  can  be  impor- 
tant for  amplitudes  as  small  as  3%. 

The  results  shown  in  Figure  36  suggest  that  a negative  value  of 
b should  be  used  in  the  approximate  analysis.  Plots  of  decay  rate  versus 
time  obtained  with  the  approximate  model  are  shown  in  Figure  37  for  linear 
combustion  (b  = 0)  and  nonlinear  combustion  (b  =-0.9).  Comparison  with 
Figure  36  shows  that  for  the  case  of  linear  combustion  driving  the  appro- 
ximate analysis  yields  larger  values  of  the  decay  rate  than  the  "exact" 
analysis.  In  order  to  assess  the  usefulness  of  the  heuristic  nonlinear 
combustion  model,  the  change  in  decay  rates  due  to  combustion  nonlinearities 
rather  than  the  actual  values  of  the  decay  r'^tes,  will  be  used  as  a basis 
for  comparing  the  approximate  and  "exact"  solutions.  For  b = -0.9  the  decay 
rate  increases  by  -4.6  sec  ^ for  the  3%  disturbance  and  by  -14.5  sec  ^ for 
the  10%  disturbance.  This  increase  in  decay  rate  is  roughly  in  proportion 
to  the  disturbance  amplitude,  as  expected  from  the  form  of  the  heuristic 
nonlinear  combustion  model  (see  Equation  (77)).  In  contrast,  the  "exact" 
analysis  yields  a -8.3  sec  ^ increase  in  the  decay  rate  for  the  3%  distur- 
bance and  a -11.7  sec  ^ increase  for  the  10%  disturbance  where  the  nonlinear 
combustion  model  is  used.  This  result  suggests  that  the  nonlinear  combus- 
tion effect  is  not  proportional  to  the  first  power  of  the  oscillation  ampli- 
tude, as  has  been  assumed  in  the  heuristic  model,  but  is  proportional  to 
some  positive  power  p of  the  amplitude,  where  p<  1. 

It  was  postulated  above  that  frequency  is  the  primary  factor  respon- 
sible for  the  differences  in  the  nonlinear  combustion  characteristics  ob- 
served in  the  two  cases  considered  in  this  study.  To  check  this  hypothesis, 
a third  case  was  also  considered  in  which  the  motor  was  shortened  from  0.597m 
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Figure  37.  Effect  of  Combustion  Nonlinearities  Upon  Decay 

Rate  for  Motor  with  2.5  p,  Particles  Using  Approx- 
imate Analysis 
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(1,958  ft)  to  0.501  m (1.642  ft)  for  the  case  with  2.5pi  particles.  This 
gave  an  oscillation  frequency  of  1079  Hz  which  is  only  slightly  higher  than 
the  frequency  obtained  in  the  longer  motor  when  particles  are  absent.  The 
corresponding  frequency  parameter  was  Q = 5.00  which  gave  a linear  response 
factor  for  the  fundamental  mode  of  “ 3.505.  In  addition  the  mean  flow 

Mach  number  at  the  nozzle  entrance  was  reduced  to  0.0654  due  to  the  decrease 
in  the  surface  area  of  the  burning  propellant.  Both  approximate  and  "exact" 
models  were  used  to  obtain  plots  of  decay  rate  versus  time  for  this  case, 
which  are  shown  in  Figure  38.  In  the  approximate  calculations,  the  same 
value  of  b was  used  as  in  the  previous  case  in  which  particles  were  absent 
(i.e.,  b = 0.32).  In  contrast  to  the  results  presented  in  Figure  36  and 
37  for  f = 918  Hz,  Figure  38  shows  a decrease  in  damping  when  the  nonlinear 
combustion  model  is  used  for  f = 1079  Hz.  Thus  combustion  nonlinearities 
increase  the  amount  of  combustion  driving  in  this  case,  which  agrees  with 
the  trend  obtained  in  the  first  case  considered  (no  particles,  1071  Hz). 

On  the  basis  of  the  "exact"  solutions,  it  appears  that  a somewhat  larger 
value  of  b should  be  used  in  the  approximate  analysis;  this  may  be  due  to 
the  slightly  higher  frequency  obtained  in  this  case  (1079  Hz  as  compared 
to  1071  Hz  obtained  previously).  These  results  indicate  that  the  effects 
of  combustion  nonlinearities  upon  motor  stability  are  strongly  dependent 
upon  the  frequency  of  oscillation. 

It  has  been  shown  in  this  section  that  for  some  cases  combustion 
nonlinearities  increase  the  amount  of  combustion  driving  as  the  amplitude 
increases.  In  the  first  case  considered  (f  = 1071  Hz,  no  particles),  the 
motor  was  unstable  to  small  amplitude  disturbances,  and  the  burning  rate 
nonlinearities  led  to  a larger  limiting  amplitude  of  the  resulting  oscil- 
lations. In  the  third  case  (f  = 1079  Hz,  2 . 5 particles)  the  motor  was 
stable  to  moderate  amplitude  disturbances,  and  the  increased  combustion 
driving  due  to  burning  rate  nonlinearities  resulted  in  smaller  decay  rates. 
The  latter  case  suggests  the  possibility  that,  if  the  combustion  nonlinear- 
ities are  sufficiently  strong,  larger  amplitude  disturbances  will  grow  in 
a motor  which  is  stable  for  small  amplitude  oscillations.  Such  a motor 
is  said  to  be  subject  to  pulsed  or  triggered  instabilities. 

The  possibility  of  pulsed  instability  was  investigated  with  the 
approximate  analysis  using  the  heuristic  nonlinear  combustion  model.  For 
the  second  case  considered  above  (918  Hz,  2.5  p.  particles),  a hypothetical 
combustion  response  was  considered  by  assuming  that  b = 3.0  (rather  than 
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the  negative  value  indicated  by  the  "exact"  model).  Solutions  were  then 
generated  for  several  different  initial  disturbance  amplitudes  using  both 
linear  and  nonlinear  burning  rate  models.  The  results  of  these  computations 
are  shown  in  Figure  39.  For  the  1%  and  3%  disturbances,  combustion  non- 
linearities  yield  smaller  decay  rates  as  expected,  while  the  10%  and  20% 
disturbances  grow  when  the  nonlinear  combustion  model  is  used.  For  the 
20%  disturbance  it  appears  that  the  oscillations  are  approaching  a limiting 
amplitude  of  about  30%.  This  is  clearly  an  example  of  pulsed  instability 

in  a linearly  stable  motor.  ; 

To  further  explore  this  case,  a plot  of  growth/decay  rate,  ^ , 
versus  amplitude  was  constructed  for  these  solutions.  This  graph  is  presented 
in  Figure  AO  which  shows  that  a increases  monotonical ly  from  negative 
values  (decay)  at  small  amplitudes,  vanishes  at  about  4%  amplitude,  and 
reaches  maximum  positive  values  (growth)  at  about  20%  amplitudes.  For 

amplitudes  larger  than  20%,  a falls  rapidly  to  zero.  The  vanishing  of  o at  i 

4%  amplitude  is  very  significant,  since  it  indicates  a threshold  amplitude 
below  which  all  oscillations  decay  and  above  which  all  oscillations  grow; 
this  amplitude  is  referred  to  as  an  unstable  limit-cycle  or  triggering 
limit.  The  rapid  drop  in  o for  large  amplitude  oscillations  indicates 
approach  of  the  solution  to  a final  limiting  amplitude  or  stable  limit  cycle 
at  which  a = 0.  At  this  point  the  increased  combustion  driving  at  higher 
amplitudes  is  balanced  by  the  increased  losses  due  to  gasdynamic  mode-coupl- 
ing which  limits  the  final  amplitude  attained.  Although  this  case  is  purely 
hypothetical,  it  illustrates  the  essential  features  of  pulsed  instability 
in  actual  rocket  motors  and  demonstrates  the  feasibility  of  using  the  appro- 
ximate nonlinear  combustion  model  to  investigate  this  phenomenon. 

In  summary,  this  study  has  shown  that  pressure  coupled  combustion 
non  1 inear  it ies  are  important  in  many  cases;  the  "exact"  calculations  indi- 
cate that  they  should  be  included  whenever  the  oscillation  amplitudes  exceed 
the  7 - 10%  range  (and  sometimes  for  amplitudes  as  small  as  3%).  Depending 
on  the  frequency  of  oscillation  and  the  various  propellant  properties,  com- 
bustion nonlinearities  may  increase  or  decrease  the  combustion  driving  as  j 

the  wave  amplitude  increases  thereby  affecting  the  growth/decay  rates  and 

limiting  amplitudes.  It  has  been  shown  that  the  heuristic  nonlinear  com-  I 

bustion  model  used  in  the  approximate  analysis  accounts  for  many  of  the  , 

effects  of  nonlinear  combustion  driving,  but  additional  information  is  needed 
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Figure  40,  Decay  Rate  vs  Amplitude  for  Nonlinear  Combustion  with 


to  determine  the  nonlinear  combustion  parameter  b.  Finally  it  has  been 
shovm  that  burning  rate  nonlinearities  may  be  an  important  factor  contribut- 
ing to  the  development  of  pulsed  instabilities,  a phenomenon  which  is  not 
predicted  on  the  basis  of  gasdynamic  nonlinearities  (mode-coupling)  alone. 

The  treatment  of  combustion  nonlinearities  presented  in  this  section 
is  by  no  means  exhaustive;  much  work  remains  to  be  done  which  is  beyond 
the  scope  of  this  project.  In  the  first  place,  the  effects  of  combustion 
nonlinearities  should  be  investigated  over  a greater  range  of  parameters 
(i.e.,  frequency,  mean  pressure,  and  propellant  properties)  than  was  possi- 
ble here.  The  heuristic  nonlinear  combustion  model  should  be  improved  to 
account  for  the  nonlinear  dependence  of  combustion  driving  upon  amplitude 
noted  earlier.  Finally,  the  most  difficult,  and  probably  the  most  impor- 
tant, work  needed  in  this  area  is  to  incorporate  the  effects  of  velocity- 
coupled  combustion  nonlinearities  into  the  "exact"  and  approximate  models. 

4 . 4 Effect  of  Nonlinear  Particle  Damping 

In  the  two  preceding  sections,  the  effects  of  nonlinearities  in 
the  gasdynamic  mode-coupling  and  the  pressure-coupled  combustion  response 
were  investigated.  In  this  section  the  effects  of  nonlinearities  in  the 
viscous  interaction  between  the  particles  and  gas  are  considered.  Using 
the  Kooker-Zinn  "exact"  model  with  the  nonlinear  particle  drag  law  (Equation 
(99)), solutions  were  calculated  for  several  particle  sizes  for  cases  with 
and  without  combustion  driving.  The  nonlinear  particle  damping  solutions 
were  then  compared  with  the  corresponding  solutions  obtained  with  the  linear 
Stokes  drag  law  (Equation  (7))  to  assess  the  importance  of  particle  damping 
nonlinearities.  In  these  studies,  nonlinear  gasdynamic  mode-coupling  was 
also  included  in  the  analysis,  but  the  pressure-coupled  combustion  response 
was  assumed  to  be  linear.  Approximate  solutions  were  also  obtained  using 
the  heuristic  nonlinear  particle  damping  law  given  by  Equation  (81),  and 
these  solutions  were  compared  with  the  corresponding  "exact"  solutions  in 
order  to  assess  the  validity  of  the  heuristic  approach. 

"Exact"  Solutions.  The  first  case  considered  was  a hypothetical 
motor  for  which  the  propellant  is  insensitive  to  the  pressure  oscillations 
(i.e.,  W = 0),  but  steady-state  combustion  is  present.  This  case  was 
chosen  in  order  to  retain  the  mean  flow  contribution  to  the  nonlinear  par- 
ticle damping,  an  effect  which  is  not  present  in  the  case  of  a particle/gas 
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mixture  in  a closed-ended  box.  Using  the  "exact"  computer  program  with 
the  nonlinear  particle  drag  option,  solutions  were  calculated  for  2 . 5p,  S o 
S 40u  for  C = 0.20,  M = 0.078,  and  f = 1071  Hz.  Plots  of  decay  rate 
versus  particle  diameter  for  nonlinear  drag  (2.5t  and  5%  amplitude)  are 
compared  with  the  corresponding  curve  for  linear  drag  in  Figure  41. 

Figure  41  shows  that  the  dependence  of  decay  rate  upon  particle 
size  for  the  nonlinear  drag  calculations  is  similar  to  that  obtained  with 
the  linear  drag  law;  both  exhibit  an  optimum  particle  size  for  maximum 
damping.  However,  the  particle  drag  nonlinearities  shift  the  peaks  to 
larger  particle  sizes.  For  the  nonlinear  drag  law,  the  maximum  damping 
occurs  at  a = 13  |j,  for  2.5%  amplitude  and  at  o = 16  p,  for  5.0%  amplitude 

as  compared  to  a ~ 8 . 5 p,  for  the  Stokes  drag  law.  The  nonlinear  particle 

drag  effect  increases  as  the  amplitude  of  the  oscillation  increases,  and 
it  also  increases  with  increasing  particle  size.  Both  of  these  effects 
are  a consequence  of  the  higher  Reynolds  numbers  resulting  from  the  in- 
crease in  particle  size  and  relative  velocity.  For  a < 10 pL  the  decay  rate 
decreases  with  increasing  amplitude,  while  damping  increases  with  increasing 
amplitude  for  a > 10  pi  . For  g < 3 pi  the  effect  of  particle  nonlinearities 
is  small  and  can  be  neglected,  but  the  nonlinear  drag  effects  become  very 
important  for  G > 15  pi  • 

The  next  case  considered  was  the  basic  set  of  motor  parameters  used 

•k 

in  the  mode-coupling  study:  = -137.0  cal/gm  (A  = 6.0,  B = 0.58,  n = 

0 575),  M = 0.078,  and  f = 1071  Hz.  Calculations  were  made  for  a = 2 . 5 p,  , 

e g 

8.0  p.  , and  20.0  p.  for  a particle  loading  of  = 0.36  using  both  linear 
and  nonlinear  particle  drag  models.  The  results  of  these  calculations  are 
shown  in  Figures  42,  43,  and  44. 

For  2.5p  particles  the  growth  to  limiting  amplitude  of  a 10%  initial 
disturbance  is  shown  in  Figure  42  for  both  linear  and  nonlinear  drag  laws. 
Particle  nonlinearities  increase  the  limiting  amplitude  from  about  9.2% 
to  about  10.5%,  a modest  effect.  For  small  particles  this  result  is  expected 
from  the  small  decrease  in  particle  damping  due  to  nonlinear  drag  effects 
at  higher  amplitudes  (Figure  41). 

For  8.0  p,  particles  the  motor  is  stable  due  to  the  greatly  increased 
particle  damping.  To  determine  the  effect  of  particle  drag  nonlinearities 
upon  the  decaying  oscillations,  decay  rate  was  plotted  as  a function  of 
amplitude  (zero-to-peak)  for  each  cycle  during  the  decay.  These  plots  for 
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Figure  42,  Effect  of  Particle  Drag  Nonlinearities  Upon 
Limiting  Amplitude  for  Motor  With  2.5  u. 
Particles  by  "Exact"  Analysis 


10%  and  20%  initial  amplitudes  for  nonlinear  drag  and  for  10%  initial  ampli- 
tude for  linear  drag  are  shown  in  Figure  43.  For  nonlinear  drag,  the  decay 
rate  increases  as  the  oscillation  decays,  approaching  a constant  value  of 
about  -273  sec  ^ for  small  amplitudes.  This  limiting  value  is  significantly 
smaller  than  the  nearly  constant  value  of  -323  sec  ^ obtained  with  the  linear 
drag  law.  There  is  very  good  overlap  between  the  nonlinear  solutions  for 
10%  and  20%  initial  amplitudes,  indicating  only  a slight  dependence  of  decay 
rate  on  the  history  of  the  oscillation. 

Similar  plots  are  shown  in  Figure  44  for  20m.  particles,  which  also 
yields  a stable  motor.  Here  the  particle  drag  nonlinearities  have  an  effect 
on  the  decay  rate  opposite  to  that  observed  for  the  8m  particles.  The 
decay  rate  decreases  as  the  wave  decays  and  the  limiting  value  of  the  decay 
rate  at  small  amplitude  (-293  sec  is  considerably  larger  than  the  value 
(-137  sec  obtained  with  the  linear  drag  law.  In  addition  there  is  a 
significant  dependence  of  decay  rate  upon  the  prior  history  of  the  oscilla- 
tion, indicated  by  the  differences  between  the  10%  and  20%  initial  ampli- 
tude solutions. 

The  nonlinear  behavior  described  above  can  be  explained  on  the  basis 
of  the  variation  of  linear  particle  damping  with  particle  diameter  a and 
the  fact  that  the  nonlinear  contribution  to  K is  always  positive  (equation 
(100)).  Thus  for  all  particle  sizes  the  effect  of  including  the  nonlinear 
drag  term  is  similar  to  the  effect  obtained  by  increasing  K without  including 
the  nonlinear  drag  term.  For  small  particles  (large  K)  increasing  K de- 
creases the  damping,  while  for  large  particles  (small  K)  increasing  K in- 
creases the  damping.  For  this  reason  including  the  nonlinear  drag  term 
decreases  the  decay  rate  for  small  particles  and  increases  the  decay  rate 

for  large  particles.  The  increasing  importance  of  particle  drag  nonlinearities 

-2  2/3 

with  increasing  particle  size  arises  because  K “ ct  and  « a 

Thus  for  sufficiently  large  a,  the  nonlinear  contribution  becomes 

larger  than  the  linear  contribution  K (Table  4). 

The  principal  factor  influencing  the  nonlinear  drag  term  is  the  abso- 
lute value  of  the  relative  velocity  between  particles  and  gas.  The  relative 
speed  is  composed  of  a steady-state  part  and  a fluctuating  part.  Due  to 
the  steady-state  lag  of  the  particles  relative  to  the  gas  in  a motor  with 
burning  at  the  lateral  boundary  (Equation  (41)),  the  steady-state  part  is 
significant.  The  steady-state  contribution  to  the  relative  velocity  ex- 
plains the  fact  that  the  damping  at  vanishingly  small  amplitude  with  non- 


116 


Zero-to-Peak  Pressure  Amplitude  (%) 

Figure  43.  Effect  of  Particle  Drag  Nonlinearities  Upon  Decay  Rates  for 


Decay  Rate,  a ( sec 


Table  4.  Linear  and  Nonlinear  Drag  Constants  Versus  o 


a 

K 

Sl 

2.5 

102.5 

0.742 

3.0 

71.18 

0.838 

4.0 

40.04 

1.015 

6.0 

17.80 

1.330 

8.0 

10.01 

1.611 

10.0 

6.406 

1.869 

15.0 

2.847 

2.449 

20.0 

1.602 

2.967 

25.0 

1.025 

3.443 

30.0 

0.712 

3.888 

40.0 

0.400 

4.710 

linear  drag  differs  from  the  corresponding  damping  with  linear  drag.  Thus 
nonlinear  particle  drag  is  important  even  for  infinitesimally  small  ampli- 
tude oscillations.  The  fluctuating  contribution  to  the  nonlinear  drag 
accounts  for  most  of  the  variation  of  decay  rate  with  amplitude  when  the 
drag  nonlinearities  are  taken  into  account. 

To  determine  the  effect  of  particle  drag  nonlinearities  upon  the 

limit-cycle  amplitudes  for  motors  with  particles  much  larger  than  2.5 p , it 

was  necessary  to  increase  the  combustion  response  and  decrease  the  particle 

* 

concentration  so  that  instability  could  occur.  Thus  Q was  increased  to 

s 

-139.4  cal/gm  giving  A = 6.00,  B = 0.55  and  n = 0.552,  and  was  decreased 
to  0.20.  Results  calculated  with  the  "exact"  analysis  for  3.0  p,  and  20 p 
particles  are  shown  in  Figures  45,  46,  and  47. 

Figure  45  shows  growth  rate  versus  number  of  cycles  for  37.  and  107. 
initial  disturbances  for  a motor  with  8.0 p particles.  This  case  is  inter- 
esting because  it  shows  that  particle  drag  nonlinearities,  like  combustion 
nonlinearities,  may  lead  to  pulsed  instabilities.  In  this  example  3% 
initial  disturbances  decay  at  a rate  of  roughly  -25  sec  while  107.  initial 
disturbances  grow  at  a rate  of  about  15  sec  This  behavior  is  a conse- 

quence of  the  nonlinear  drag  law  which  causes  a decrease  in  particle  damping 
with  increasing  amplitude  (Figure  43)  allowing  the  combustion  process  to 
drive  oscillations  for  sufficiently  large  amplitudes. 
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The  particle  nonl  inear  it  ies  have  the  opposite  effect  for  20  (i,  par- 
ticles as  shown  in  Figure  46.  Using  the  linear  drag  law,  the  analysis 
predicts  growth  of  oscillations  to  a limiting  amplitude  exceeding  8%.  For 
the  nonlinear  drag  law,  the  greatly  increased  particle  damping  causes  the 
oscillations  to  decay  to  zero  (or  at  least  to  a limiting  amplitude  less 
than  1.5%).  The  corresponding  pressure  waveforms  for  this  case  are  pre- 
sented in  Figure  47,  which  shows  the  dramatic  effect  of  particle  drag  non- 
linearities  upon  the  amplitude  of  the  oscillations  but  little  or  no  effect 
on  the  shape  of  the  pressure  waveform. 

Approximate  Solutions.  In  order  to  assess  the  validity  of  the 
heuristic  nonlinear  particle  drag  law  described  in  Section  2.8,  approximate 
solutions  were  obtained  for  8 n and  20  p,  particles  and  were  compared  with 
the  corresponding  "exact"  solutions.  These  results  are  presented  as  plots 
of  decay  rate  versus  amplitude  in  Figures  48  and  49  which  correspond  to  the 
"exact"  solution  presented  in  Figures  43  and  44. 

Figure  48  shows  approximate  solutions  for  decay  rate  as  a function 
of  amplitude  for  a motor  with  8 ij,  particles.  For  nonlinear  drag  these 
curves  are  shown  for  several  values  of  the  parameter  C which  appears  in 
Equation  (81).  For  the  case  of  C = 0 the  amplitude-dependent  part  of  the 
nonlinear  drag  vanishes  and  only  the  steady-state  contribution  remains. 

For  positive  values  of  C the  decay  rate  increases  as  the  amplitude  of  the 
wave  decreases,  which  agrees  with  the  "exact"  calculations.  In  addition 
the  approximate  solutions  agree  with  the  "exact"  solutions  regarding  the 
prediction  that  the  particle  damping  at  low  amplitudes  is  smaller  when  the 
nonlinear  particle  drag  law  is  used  (due  to  the  steady-state  contribution 
to  the  Reynolds  number). 

In  comparing  the  approximate  and  "exact"  solutions  shown  in  Figures 
43  and  48,  it  should  be  noted  that  the  two  analyses  predict  different  decay 
rates  for  linear  particle  drag  (i.e.,  o = -323  sec”^  lor  the  "exact"  analysis 
compared  to  a = -400  sec  ^ for  the  approximate  analysis).  Thus  quanti- 
tative assessment  of  the  validity  of  the  postulated  nonlinear  drag  model 
must  be  made  by  comparing  the  "exact"  and  approximate  results  on  the  basis 
of  the  change  in  decay  rate  due  to  the  addition  of  particle  drag  nonlinear- 
ities rather  than  on  the  basis  of  the  actual  values  of  the  decay  rates. 
Therefore  two  quantities,  Aa  and  Aa  2 considered  in  comparing 
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tho  approxiniatt-  and  "exact"  solutions.  Ttie  first  quantity,  La  ^ is  the 
difference  between  the  decay  rate  at  small  amplitude  computed  using  the 
nonlinear  drag  law  and  the  corresponding  decay  rate  obtained  using  the 
linear  drag  law;  thus,  ~ ainonlinear  drag,  small  amplitude)  - a 

(linear  drag).  By  comparing  the  approximate  values  of  with  the  "exact" 

values,  the  validity  of  the  steady-state  term,  Kgg . in  Equation  (81)  can 
be  assessed.  The  second  quantity,  2'  difference  between  the  decay 

rate  at  finite  amplitude  and  the  decay  rate  at  vanishingly  small  amplitude 
as  calculated  using  the  nonlinear  drag  models;  thus,  “ a(finite 

amplitude)  - a(small  amplitude).  Comparing  the  approximate  and  "exact" 
values  of  La 2 for  a given  amplitude  can  be  used  to  determine  the  proper 
value  of  the  parameter  C appearing  in  the  unsteady  contribution  to  the  non- 
linear particle  drag  (Equation  (81)). 

For  the  case  of  8.0^  particles  shown  in  Figure  48,  values  of  a^>j 

and  AO,  were  calculated  using  both  approximate  and  "exact"  analyses.  For 

^ . . -1  -1 
Aetj^  the  approximate  analysis  yields  35  sec  compared  to  50  sec  for  tlie 

"exact"  analysis.  Thus  the  heuristic  nonlinear  drag  model  underestimates 

the  steady-state  contribution  to  the  nonlinear  drag  for  8 particles. 

At  4%  amplitude  the  approximate  analysis  yields  Ao,  = 87  sec  ^ for  C = 1 

and  AO.,  = 48  sec  for  C = 0.5,  while  the  "exact"  analysis  yields  ao, 

^-1  . . ' " 

= 63  sec  . This  result  indicates  that  a value  of  C = 0.7  would  yield 

better  agreement  with  the  "exact"  analysis  for  the  case  of  8,1  particl-  s. 

Figure  49  shows  similar  decay  rate  versus  amplitude  curves  obtained 

with  the  approximate  analysis  for  a motor  with  20 j.  particles.  Again  the 

qualitative  agreement  with  the  "e.tact"  solutions  (Figure  44)  is  good;  the 

approximate  solutions  (0  = 1),  also  exhibit  a decreasing  decay  rate  as  the 

amplitude  increases  and  a limiting  decay  rate  at  low  amplitude  wiiich  is 

considerably  larger  than  the  value  obtained  with  the  linear  drag  law.  For 

this  case  the  approximate  method  gave  La.  = -205  sec  ^ as  compared  to 
- 1 

A®  = -156  sec  from  the  "exact"  analysis.  Thus  the  approximate  analysis 
overestimates  the  steady-state  contribution  to  the  nonlinear  drag  for  20  n, 
particles,  which  is  opposite  from  the  discrepancy  obtained  for  8 p,  particles. 
Also  the  amplitude  dependent  part  of  the  nonlinear  drag  is  overestimated 

-1 

for  C = 1;  at  4%  amplitude  the  approximate  analysis  yielded  l<^2  ~ “160  sec 
while  the  "exact"  analysis  (10%  initial  amplitude)  gave  AQ'2  = -117  sec”^. 


125 


As  in  the  case  of  8^  particles,  these  results  indicate  that  C = 0.7  will 
give  better  agreement  regarding  the  amplitude  dependence  of  the  nonlinear 
particle  damping. 

These  results  presented  above  indicate  that  the  heuristic  nonlinear 
drag  model  used  in  the  approximate  analysis  gives  a reasonable  qualitative 
description  of  the  effects  of  particle  drag  nonlinearities  on  the  decay 
rates  for  stable  solid  rocket  motors.  To  obtain  better  quantitative  agree- 
ment between  the  approximate  and  "exact"  nonlinear  particle  damping  so- 
lutions, a more  extensive  parametric  study  is  needed  to  determine  how  the 
heuristic  nonlinear  particle  drag  model  can  be  improved.  Such  a study 
was  beyond  the  scope  of  this  project  and  is  recommended  for  future 
research. 

4 . 5 Solutions  for  T-Burners 

In  this  subsection  the  Galerkin  method  is  used  to  obtain  approxi- 
mate solutions  for  a few  selected  T-burner  configurations  in  order  to 
demonstrate  the  applicability  of  the  approximate  analysis  to  T-burners. 

For  most  of  the  cases  considered  here,  the  length,  bore,  propellant  pro- 
perties, burning  rate  and  combustion-product  properties  (gaseous  and 
particulate)  for  the  T-burner  are  assumed  to  be  the  same  as  those  for  the 
laboratory  pulse  motor  considered  in  Section  4.2.  Thus  the  dimensionless 
radius  for  the  T-burner  is  R = 0.051  and  the  dimensionless  velocity  of 
the  combustion  products  leaving  the  burning  propellant  surfaces  is 
Uj^  = 0.002.  The  transient  combustion  parameters  are  the  same  as  those 
used  for  the  motor  solutions  with  particles;  that  is,  A = 6.00,  B = 0.58, 
and  n = 0.575.  The  dimensionless  frequency  parameter  is  assumed  to  be 
fi  = 4.2  which  yields  a response  factor  of  ” 4.14.  The  dimensionless 

length  of  the  center  vent  is  given  by  8^  = 0.1,  while  the  dimensionless 
plug  *^low  length  for  the  pipe  connecting  the  burner  with  the  surge  tank 

is  assumed  the  value  L = 0.166.  In  accordance  with  the  theoretical 

ef  f 

limitations  discussed  in  Section  2.6,  flush  grains  (S  /S  = 1.00)  are 

assumed  for  all  cases  considered  in  this  study. 

In  obtaining  the  approximate  T-burner  solutions  the  following 

parameters  were  varie<1:  (1)  the  length  of  the  cylindrical  grains, (2)  the 

combustion  response  f , (3)  the  particle  concentration  C , and  (4)  the 

1 m 

vent  factor  . Grain  configurations  considered  were  end-burning  only 
(8=0)  and  cup  grains  < 3 = 0.1  and  0.3).  The  effect  of  combustion 
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driving  on  the  T-burner  solutions  was  determined  by  comparing  the  growth 
rate  for  R ^ =4.14  with  the  decay  rate  obtained  when  the  burning  rate 

is  insensitive  to  pressure  oscillations  ( R =0).  Solutions  for  the  case 
of  no  particles  (C^  = 0)  were  also  compared  with  those  obtained  for  the 
case  of  2.5u  -diameter  particles  with  C = 0.36.  The  effect  of  the  center 
vent  was  assessed  by  considering  the  case  of  vent  gain  (V.  = 0)  and  no 
vent  effect  (V^  = 1). 

The  case  of  end-burning  (disk  grains,  g = 0)  was  considered  first. 
Figure  50  shows  growth/decay  rates  for  3%  initial  IL-mode  disturbances 
in  end-burning  T-burners  for  four  cases;  (1)  no  particles,  R = 0,  (2) 

no  particles,  R^  = 4.14,  (3)  2.5p,  particles,  R = 0,  and  (4)  2.5p,  par- 
ticles, R^  = 4.14.  In  each  case  both  linear  (IL  mode  only)  and  nonlinear 
(5  modes)  solutions  are  presented.  For  the  first  case,  the  nonlinear 
solutions  decay  at  approximately  the  linear  rate  of  -1.8  sec  ^ during  the 
first  few  cycles;  the  slight  amount  of  damping  is  due  to  the  effect  of 
the  mean  flow  and  the  center  vent.  Similarly  the  initial  growth  rate  of 
the  nonlinear  solutions  for  the  second  case  (with  combustion  driving)  is 
close  to  the  linear  value  of  41  sec  In  both  cases  the  rapid  rise  in 

growth  rate  after  seven  cycles  is  a result  of  strong  waveform  steepening 
due  to  nonlinear  gasdynamical  mode-coupling  (generation  of  higher  harm- 
onics). When  2.5m.  particles  are  present,  the  T-burner  pressure  oscil- 
lations decay  at  a nearly  constant  rate  of  -68  sec  ^ for  R = 0 (case  3) 
and  -37  sec  ^ for  = 4.14  (case  4)  due  to  the  increased  attenuation 

of  the  higher  harmonics  by  the  particles. 

In  order  to  obtain  growth  of  oscillations  to  limiting  amplitude 
for  T-burners  with  metallized  propellants  it  is  often  necessary  to  increase 
the  surface  area  of  the  burning  propellant  by  using  cup  grains.  This  is 
the  situation  in  the  case  selected  for  this  study,  for  a cup-length  of 
8 = 0.1  is  necessary  to  obtain  spontaneous  growth  of  oscillations  when 

the  2.5  ^ particles  are  present.  Growth  rates  for  cup-lengths  of  8 = 0 
(end-burning  only),  8 “ 0.1,  and  8 “ 0-3  are  compared  in  Figure  51.  For 
8 =0.1,  the  growth  rate  of  the  nonlinear  solutions  rapidly  decreases 

from  the  linear  value  of  12  sec  ^ to  nearly  zero,  indicating  approach  to 
limiting  amplitude  (about  7%  in  this  case).  For  8 * 0.3,  the  initial 
growth  rate  is  much  larger  (about  100  sec  ^)  than  for  8 “ 0.1,  but  it  de- 
creases rapidly  as  the  amplitude  increases.  For  the  latter  case,  the 
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solution  was  not  carried  out  to  limiting  amplitude,  but  it  probably  exceeds 

20%. 


The  importance  of  the  center  vent  for  this  case  was  determined 
by  comparing  solutions  for  = 0 (vent  gain)  with  the  corresponding  solu- 

tions for  = 1 (no  vent  effect).  These  comparisons  were  made  for  both 
end-burning  and  cup-grain  configurations,  for  T-burners  without  combustion 
driving  ( rt  = 0)  or  particle  damping.  The  results  of  these  calculations 
are  given  in  Table  5.  The  quantity  ^ is  the  difference  between  the 
decay  rates  obtained  with  =0  and  = 1;  thus,  represents  the 

vent  gain  due  to  the  raean-f low/acoustic  interaction.  From  Table  5 it  is 
seen  that  AQ'.^  is  a rather  small  quantity  that  increases  as  g increases 
due  to  the  resulting  increase  in  steady-state  flow  velocity  in  the  vent 
region.  These  results  indicate,  that,  at  least  for  the  cases  considered, 
the  raean-f low/acoust ic  vent  effect  is  relatively  unimportant  and  may  be 
neglected  in  the  nonlinear  analysis  of  T-burners. 


Table  5,  Effect  of  Center  Vent 


— 

Case 

— 

3 

a,  (sec‘^) 

V,  = i 

1 

0.0 

-1.83 

-1.76 

2 

0.1 

-22.18 

-21.98 

3 



0.3 



-59.03 



-58.55 

0.07 

0.20 

0.48 


Pressure  waveforms  and  mode-ampl  i tudt*  functions  wrf  obtained  for 
two  of  the  T-burner  cases  considered  in  this  study.  Ihesi-  results  are 
shown  in  Figure  52  for  an  end-burning  T-burner  without  part  . ,les  and  in 
Figure  53  for  a cup  grain  configuration  (r  = d.l'  wili.  2 d araeter 

particles;  in  both  cases  a five-mode  senes  was  used  and  • .esponse 
function  was  = 4.1A. 

In  the  first  case  (Figure  52),  a lOJ  amplitud>  1 , -mode  initial 
disturbance  rapidly  grows  and  steepens  after  two  wave  les  as  higher 
harmonics  are  generated  by  nonlinear  gasdynamic  roode-t oupl i ng . As  in  the 
case  of  solid  rocket  motors  (Figure  8),  the  amplitude  of  the  fundamental 
mode  is  considerably  larger  than  the  amplitudes  of  the  higher  harmonics 
(the  response  functions  for  the  higher  modes  are  considerably  smaller  than 
the  response  function  for  the  fundamental  mode).  It  is  also  seen  that 
the  amplitudes  of  the  higher  modes  decrease  as  the  axial  mode  number  in- 
creases, except  that  the  5L  mode  has  a larger  amplitude  than  the  4L  mode. 
The  latter  result  is  probably  due  to  the  finite  number  of  modes  (i.e., 
five)  used  in  the  series,  and  it  accounts  for  the  fifth  harmonic  "wiggles" 
apparent  in  the  pressure  waveform. 

In  the  second  case  (Figure  53)  when  particles  are  present,  a 7.5% 
amplitude  1-L  mode  initial  disturbance  is  shown  developing  into  a nearly 
sinusoidal  limiting-amplitude  pressure  waveform.  In  contrast  to  the  pre- 
vious case  the  amplitudes  of  the  higher  harmonics  decrease  rapidly  as  the 
mode  number  increases;  in  fact  it  appears  that  the  3L,  4L  and  5L  modes 
are  negligible  and  that  a two  or  three  mode  series  is  adequate  to  represent 
the  solution. 

In  order  to  adequately  evaluate  the  applicability  of  the  approxi- 
mate technique  to  the  nonlinear  analysis  of  T-burners,  the  approximate 
T-buroer  solutions  were  compared  with  available  "exact"  solutions  as  well 
as  with  experimental  data.  Since  time  did  not  permit  extension  of  the 
Kooker-Zinn  "exact"  analysis  to  analyze  T-burner  configurations,  Levine's 
T-burner  calculations^  were  used  in  these  comparisons.  The  case  selected 
for  comparison  with  the  "exact"  analysis  was  for  = 7.06  (Table  A-1 

of  Reference  5)  for  which  0 = 0.216  (tubular  grains  without  end-bur. .ing 
disks).  The  combustion  parameters  used  by  Levine  were  A = 8.8,  B = 0.67, 
n = 0.3,  and  0 = 7.80  which  gave  a response  function  of  if  ^ = 2.03. 

The  particles  were  characterized  by  an  average  diameter  of  3.0m.  and  a 
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of  Mode  Amplitude  Function,  A.(t)  Pressure  Pe 
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concentration  of  C = 0.36. 

ra 

T-burner  pressure-t ime  histories  obtained  with  the  Galerkin  method 
and  Levine's  "exact"  analysis  are  compared  in  Figure  54  for  3%  initial 
disturbances.  The  resulting  solutions  are  seen  to  disagree  regarding  the 
initial  growth  rate  and  limiting  amplitude.  The  approximate  solutions 
grow  very  slowly  as  they  appear  to  be  near  limiting  amplitude  (probably 

about  4%),  while  the  "exact"  solutions  are  still  growing  rapidly  after  ' 

12  wave  cycles  (eventually  reaching  37%  limiting  amplitude).  Both  solu- 
tions agree,  however,  regarding  the  nearly  sinusoidal  waveshape  and  the 
double  frequency  small  amplitude  oscillation  at  the  center  vent.  Compari- 
sons with  experimental  data  for  T-burners  are  presented  in  the  next 
section . 

The  results  presented  in  this  subsection  indicate  that  the  Galerkin 
method  can  be  used  to  obtain  approximate  solutions  for  growth/decay  rates 
and  limiting  amplitudes  for  flush-grain  T-burner  configurations.  The 
limited  parametric  study  performed  indicates  that  the  approximate  solutions 
exhibit  the  expected  trends  regarding  the  effect  of  cup-grain  length,  com- 
bustion response  factor,  particle  concentration,  and  center  vent  upon  the 
behavior  of  the  pressure  oscillations.  It  is  difficult  to  draw  conclusions 
regarding  the  accuracy  of  the  approximate  solutions  from  the  limited  com- 
parison with  "exact"  theory  given  above,  since  Levine's  model  differs  sig- 
nificantly from  the  present  approximate  model  (i.e.,  it  accounts  for  the  | 

thermal  particle/gas  interaction  as  well  as  nonlinear  combustion  and  par-  j 

tide  drag  nonlinearities).  Furthermore,  the  limiting  amplitude  predicted 
by  the  "exact"  analysis  is  far  higher  than  that  obtained  from  the  measured 
data  (37%  as  compared  to  8%  for  the  data).^  It  is  evident  that  much  more 
work  needs  to  be  done  to  further  develop  and  improve  both  the  approximate 
and  "exact"  T-burner  models,  and  that  more  extensive  comparisons  of  the 
approximate  T-burner  solutions  with  "exact"  solutions  and  experimental 
data  are  needed. 

4 . 6 Comparisons  With  Experimental  Data 

In  the  preceding  subsections,  the  approximate  analysis  has  been 
evaluated  on  the  casis  of  comparisons  with  available  "exact"  solutions 
for  both  motors  and  T-burners.  Since  the  usefulness  of  an  analytical 
method  lies  in  its  ability  to  predict  observed  behavior,  the  approximate 
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Figure  54,  Comparison  of  Appi 
Solutions 
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solutions  were  also  compared  with  available  T-burner  and  motor  test  data. 
The  results  of  this  study  are  presented  in  this  subsection. 

Comparison  With  Motor  Test  Data.  Di.e  to  the  unavailability  of 
any  other  motor  test  data,  the  experimental  data  considered  by  Levine  and 
Culick^  was  also  used  in  this  study.  This  data  was  obtained  by  Aerojet 
several  years  ago  for  ANB  3066  propellant  using  a small  laboratory  pulse 
motor.  This  motor  has  a grain  length  of  about  0.597  m (23.5  in)  and  an 
initial  bore  of  about  50  mm  (2  in);  it  is  the  same  motor  configuration 
considered  previously  in  Sections  4.2,  4.3,  and  4.4.  The  physical  pro- 
perties of  the  ANB  3066  propellant  and  its  gaseous  and  particulate  com- 


bustion products  have  already  been  given  in  Table  1 of  Section  4.2.  The 
motor  parameters  for  the  four  cases  considered  by  Levine  and  Culick^  are 


given  in  Table  6, 


Table  6.  Motor  Parameters  for  Laboratory  Pulse  Motor  I 


a 

] 

1 

I 

i 


The  Mach  number  at  the  nozzle  entrance,  which  is  needed  in  the 
approximate  analysis,  was  determined  from  the  ratio  of  port  area  to  throat 
area  using  the  isentropic  flow  relations  (the  low  flow  coefficient  of  the 
nozzle  was  already  taken  into  account  by  Levine  and  Culick  by  adjusting 
the  throat  area).  The  flame  temperature  T^  was  taken  from  Figure  7-17 


Nominal  Throat 
Diameter  (cm) 

Pulse 

Number 

Port  Area 
(cm^) 

Throat  Area 
(cm^) 

P 

(knt/m^) 

M 

e 

^f 

(K) 

1.97 

1 

21.5 

2.83 

10810 

0.0782 

1.97 

3 

36.3 

2.95 

14800 

0.0481 

2.20 

1 

24.4 

3.57 

8330 

0.0869 

2.20 

2 

30.5 

L-  

3.63 

9740 

0.0703 

3 


of  Reference  5.  For  2 . 0 tx  particles  the  drag  constant  K was  about  46.7 
for  all  four  cases. 

In  order  to  obtain  approximate  solutions  for  the  test  conditions 

given  above,  the  transient  burning  rate  parameters  (i.e.,  A and  B)  must 

be  estimated.  Levine  and  Culick^  made  such  estimates  based  on  the 

following  criteria:  (1)  the  response  function  must  be  greater  than  3 

for  instability,  (2)  the  values  of  A and  B must  not  result  in  nonlinear 

★ * 

intrinsic  instability,  and  (3)  the  corresponding  values  of  and  should 

be  reasonable.  For  Cases  1 and  4 this  procedure  was  relatively 

straightforward,  bit  for  Cases  2 and  3 the  considerable  spread  in  mean 

pressure  caused  difficulties.  For  the  latter  two  cases  Levine  and  Culick 

★ 

postulated  a pressure  dependence  of  above  and  beyond  that  which  is 
accounted  for  by  the  Levine-Cul ick  nonlinear  transient  burning  rate  model. 

In  addition  to  A and  B,  the  dimensionless  frequency  parameter  was 
calculated  using  the  steady-state  burning  rate  given  by 

F = 0.813  (P/3448)°’^  (108) 

where  F is  the  surface  regression  rate  in  units  of  cra/sec  and  P is  the 

2 

steady-state  pressure  in  knt/m  . The  transient  burning  rate  parameters 
assumed  by  Levine  and  Culick  were  used  as  the  starting  point  in  this  inves- 
tigation; they  are  given  in  Table  7. 


Table  7.  Transient  Burn  Rate  Parameters 


Case 

A 

B 

w 

1 

6.00 

0.55 

m 

• 

2 

5.98 

0.53 

3.65 

3 

6.02 

0.54 

5.06 

4 

6.00 

0.55 

4.3 
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Approximate  solutions  were  obtaiiied  for  each  of  the  four  cases 
described  above  using  the  Galerkin  method  with  a five-mode  series.  In 
each  case,  using  Levine  and  Culick's  estimates  of  A,  B,  and  Q along  with 
n = 0.3,  the  approximate  analysis  predicted  a decaying  oscillation.  There- 
fore the  parameters  A and  B were  varied  slightly  in  order  to  increase  the 
propellant  response  to  the  level  necessary  to  sustain  oscillations.  The 
results  of  these  calculations  are  now  presented  and  discussed  for  each 
of  these  cases  in  the  following  paragraphs. 

The  results  for  Case  1 (1.97-cm  throat.  Pulse  #1)  are  summarized 
in  Table  8.  The  measured  limiting  amplitude  for  this  case  was  about  3%, 


Table  8,  Approximate  Solutions  for  Cases  1 and  4 


while  the  approximate  analysis  predicts  a limiting  amplitude  slightly 
larger  than  5%  for  B = 0.543,  For  comparison  Levine  and  Culick  obtained 
an  amplitude  of  4.2%  with  B = 0,550, 

Case  4 (2.20-cm  throat,  Pulse  #2)  is  considered  next  because  the 

mean  pressure  is  nearly  the  same  as  for  Case  1,  allowing  the  same  values 

of  A and  B to  be  used  (Table  7).  Case  4 differs  from  Case  1 due  to  the 
smaller  value  of  (0.0703  as  compared  to  0.0782  for  Case  1).  Using 
A = 6.0  and  B = 0.543  as  before,  the  approximate  analysis  yielded  a 
limiting  amplitude  of  a little  over  2X  which  is  in  excellent  agreement 

with  the  measured  value  of  slightly  over  2%.  For  this  case  Levine  and 

Culick  obtained  a 3%  oscillation  with  B = 0.55.  This  comparison  shows 
that  the  approximate  analysis  correctly  predicts  the  experimentally  ob- 
served trend  of  decreasing  limiting  amplitude  with  decreasing  steady  state 
Mach  number. 

The  approximate  solutions  for  Case  2 (1.97-cm  throat.  Pulse  #3) 
and  Case  3 (2.20-cm  throat.  Pulse  #1)  are  presented  in  Table  9.  In  each 
case,  the  first  set  of  parameters  corresponds  to  Levine's  values,  which 
yield  decaying  oscillations.  Since  the  values  of  O for  these  cases  yield 
values  of  considerably  below  the  peak  value,  variations  of  the  para- 
meter B alone  are  ineffective  in  raising  sufficiently  to  drive  in- 

stability. Since  n is  not  a strong  function  of  pressure  and  Pi  is  fairly 
well  determined  for  the  given  chamber  pressures,  the  only  parameter  left 
which  can  be  varied  is  A.  Thus  the  parameter  A was  changed  (and  B as  well 
for  Case  2)  to  yield  a response  factor  of  about  3.8.  For  Case  2 the  second 
set  of  parameters  yields  a limiting  amplitude  of  about  1.5%  in  good  agree- 
ment with  the  measured  value  of  1-3%  (modulated).  For  Case  3 the  second 
set  of  parameters  gives  an  amplitude  in  excess  of  6%  compared  to  measured 
values  of  2-2.5%,  but  a slight  decrease  in  A should  give  better  agreement. 
It  is  difficult  to  assess  the  validity  of  the  approximate  solutions  on 
the  basis  of  these  last  two  cases,  since  the  values  of  the  transient  com- 
bustion parameters  A and  B used  in  obtaining  these  solutions  are  uncertain. 
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Table  9.  Approximate  Solutions  for  Cases  2 and  3 


' ‘n 

Case 

A 

B 

Initial 

Amplitude 

Oscillation  after  12  cycles 

Amplitude  (7«) 

O’  (sec  ^) 

2 

5.98 

0.530 

57o 

2.74 

2.5 

-55.8 

2 

5.60 

0.517 

5% 

3.90 

4.5 

-5.9 

2 

5.60 

0.517 

27o 

3.90 

1.9 

-0.8 

3 

6.02 

0.540 

57o 

2.36 

1.4 

-122 

"k 

3 

6.35 

0.540 

5% 

3.81 

5.6 

19.2 

★ 

3 

6.35 

0.54 

27o 

3.81 

2.3 

21.6 

Modulated 


Comparison  With  T-Burner  Test  Data.  The  T-burner  test  data  con- 
sidered in  this  study  are  the  same  three  cases  considered  by  Levine  and 
Culick^  for  comparison  with  the  "exact"  analysis.  This  data  was  selected 
from  the  ANB  3066  (3.8-cm  bore)  data  reported  in  Reference  22  on  the  basis 
of  two  requirements;  the  measurements  were  taken  near  the  flush  grain 
conditions  and  limiting  amplitudes  were  reported.  In  these  cases,  the 
T-burner  length  was  62.2  cm  (24.5  in),  the  steady-state  pressure  was 
3448  knt/m  (500  psi),  and  the  propellant  configuration  consisted  of 


^^Beckstead,  M.W. , et  al.  "Variable  Area  T-Burner  Investigation,"  AFRPL 
TR-72-85,  December  1972. 


tubular  propellant  grains  at  each  end  of  the  burner  without  end-burning 

disks.  The  principal  parameter  of  interest  in  this  study  was  the  area 

of  the  burning  propellant  surface  which  was  specified  by  the  ratio 

•k 

s./s  and  determined  by  the  length  of  the  propellant  grains  L,  . The 

b CO 

pertinent  data  for  the  three  values  of  S./S  selected  is  given  in 

Table  10;  for  all  cases  R = 0.0306,  0^  = 0.0542,  = 1.0  (no  vent 

effect),  and  L ,_  = 0.163.  The  particle  diameter  was  assumed  to  be 
ef  f 

3 u,  and  the  particle  concentration  was  assumed  to  be  C = 0.36. 

5 • 

Following  Levine  and  Culick  , the  combustion  response  was  assumed  to  be 
given  by  A = 8.8,  B = 0.67,  and  n = 0.3. 


Table  10.  T-Burner  Parameters  for  Comparison  Study 


Case  S,  /S 

b CO 


(Hz)  (cm/sec) 


1 1.48  0.0465  772 


0.856 


6.62  0.00322 


2 5.58  0.1714  816 


0.815  . 7.71  0.00307 


3 7.06  0.2163  820 


0.813 


7.80  0.00306 


Approximate  solutions  were  obtained  for  the  above  three  cases  using 
both  single-mode  (ID  and  five-mode  series  expansions.  The  single-mode 
results  are  presented  in  Table  11;  they  correspond  to  exponential  growth 
or  decay  of  oscillations  at  small  amplitudes.  Since  the  approximate 
analysis  predicts  more  stable  behavior  than  that  indicated  by  the  measured 
data  for  B = 0.67  (Levine's  value),  approximate  solutions  for  B = 0.635 
are  also  shown  in  Table  11  which  give  better  agreement  with  the  experi- 
mental data.  The  additional  case  of  S^/S^^  = 0 (no  burning,  no  mean  flow) 
was  also  considered;  for  this  case  the  approximate  analysis  predicted  a 


decay  rate  of  -88.2  sec  ^ as  compared  to  a measured  value  of  -77  sec 
(extrapolated) . 


Table  11,  Comparison  of  Single-Mode  Solutions  with  Measured  Data 


Case 

SJS 

b CO 

B = 

0.670 

B = 0.635 

Measured 

O' 

(sec 

a 

(sec 

O' 

(sec 

1 

1.48 

2.15 

-64.9 

3.25 

-47.6 

-50 

2 

5.58 

2.06 

-4.9 

2.83 

59.9 

54 

3 

7.06 

2.03 



15.1 

2.75 

85.0 

88 



Both  five-mode  and  single-mode  solutions  for  B = 0.635  are  compared 
with  the  measured  data  in  Figure  55,  which  shows  plots  of  growth  rate 
versus  number  cf  wave  cycles  for  3%  amplitude  initial  disturbances.  For 
S^/S  = 1.A8  the  approximate  solutions  are  in  very  good  agreement  with 

b CO 

the  measured  data,  while  reasonably  good  agreement  is  obtained  for 

S, /S  = 5.58  and  S^/S  = 7.06.  The  decline  in  growth  rate  for  the  latter 

b CO  b CO 

two  cases  after  six  cycles  indicates  that  the  oscillations  in  the  T-buiner 
are  approaching  limiting  amplitude.  Figure  56  shows  the  continuation  of 
the  calculations  for  40  cycles  for  these  two  cases,  which  yields  limiting 
peak-to-peak  amplitudes  of  about  31%  for  = 5.58  and  about  43%  for 

S,  /S  = 7.06,  corresponding  to  measured  values  of  9.0%  and  15.6%  respec- 

b CO  r o 

tively . 

Three  significant  observations  regarding  the  applicability  of  the 
approximate  T-burner  analysis  can  be  made  from  the  comparisons  presented 
above.  First,  the  linear  growth  rates  predicted  by  the  approximate 


143 


Growth  Rate,  a(sec 


Number  of  Cycles 

Figure  56.  Approach  to  Limiting  Amplitude  for 
T-Bumers  by  Approximate  Analysis 
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T-burner  analysis  do  not  agree  well  with  those  obtained  by  Levine's  "exact" 
analysis  (Levine  obtained  good  agreement  with  the  measured  daf-  for 
B = 0.67).  Second,  for  a somewhat  stronger  propellant  respons 
(B  = 0.635)  the  approximate  analysis  agrees  well  with  the  measured  data 
regarding  the  effect  of  upon  the  linear  growth  rates.  Finally, 

the  approximate  analysis  predicts  limiting  amplitudes  which  are  far  higher 
than  the  measured  values  indicate,  even  though  the  linear  growth  rates 
are  in  good  agreement  (Levine's  "exact"  model  also  predicts  excessive 
limiting  amplitudes^).  The  results  of  this  limited  comparison  study  indi- 
cate that,  although  the  present  approximate  T-burner  analysis  may  be  useful 
in  predicting  trends,  further  improvement  and  development  of  the  T-burner 
model  is  needed  before  it  becomes  an  accurate  indicator  of  nonlinear  T- 
burner  behavior. 

4 . 7 Application  of  the  Method  of  Averaging 

As  discussed  in  Section  2.4,  the  Method  of  Averaging  (MOA)  can 
be  used  to  reduce  the  computation  time  required  to  obtain  approximate 
solutions  at  the  cost  of  introducing  additional  error  into  the  analysis. 

In  order  to  determine  whether  the  savings  in  computation  time  obtained 
with  the  MOA  compensates  for  the  resulting  loss  of  accuracy,  a comparison 
study  was  conducted.  In  this  investigation,  approximate  solutions  obtained 
using  the  MOA  were  compared  with  the  corresponding  solutions  obtained  with 
the  Galerkin  method  and  the  "exact"  analysis  in  order  to  determine  the 
range  of  applicability  of  the  MOA.  The  results  of  this  study  are  presented 
in  this  subsection. 

For  cases  in  which  particles  are  not  present,  the  MOA  and  the 
Galerkin  method  are  in  excellent  agreement,  as  illustrated  by  the  waveforms 
shown  in  Figure  57.  For  cases  with  particle  damping,  however,  the  MOA 
predicts  a larger  decay  rate  than  the  Galerkin  method.  For  instance, 
solutions  were  obtained  by  both  methods  for  the  case  of  the  attenuation 
of  800  Hz  waves  in  a gas-particle  mixture  enclosed  in  a box  with  a particle 

loading  of  C = 0.36  and  a particle  diameter  of  2 . 5 p,  . In  this  case,  the 

™ -1  . 
linear  decay  rate  obtained  using  the  Galerkin  method  was  -56  sec  , while 

the  linear  decay  rate  obtained  using  the  MOA  was  found  to  be  -72  sec 

Since  previously  published  results  by  Culick^^  showed  good  agree- 
ment between  the  MOA  solutions  and  Levine's  "exact"  solutions  for  cases 
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with  particle  damping,  Culick's  equations  were  also  examined  to  aid  in 
determining  the  source  of  discrepancy.  Unlike  the  present  analysis  where 
a wave  equation  for  the  velocity  potential  is  derived,  Culick  has  formu- 
lated the  problem  by  developing  an  inhomogeneous  wave  equation  for  the 
pressure  perturbation,  p' . Considering  a single  mode  expansion  consisting 
of  the  IL  'rode,  an  amplitude  function  is  defined  such  that 


p7p^  = y^(x) 


u'  = (1/yK^)  ^1  Hi 

dt  dx 


(109) 

(110) 


where  Y describes  the  mode  structure  and  = ui^/Cj  is  the  wave  number. 
To  determine  the  inhomogeneous  wave  equation  was  multiplied  by  and 
integrated  over  the  length  of  the  chamber  (this  step  is  equivalent  to  the 
Galerkin  method)  to  obtain  the  following  second  order  equation  for  : 


where 


dt 


2 -n 
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In  this  equation. 


(u'  - u') 


Ps^ 

18  p, 


and  6Up  and  6F^  are  defined  in  Reference  10  (i.e.,  Equations  (8.15)  and 
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(8.13)).  Using  the  referenced  equations  in  Equations  (111),  the  equation 


describing  the  behavior  of  becomes; 


d\ 


dt 


1 2 

+ ‘“l  \ = 
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1 +C 
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dt 
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Td  dt 


^ / dT|, 

J ' ■*  JT' 


(112) 


Since  the  right-hand-side  of  Equation  (112)  contains  a second  deri- 
vative term,  it  is  not  in  the  proper  form  for  solution  by  the  MOA.  This 
equation  can  be  simplified  and  rewritten  in  the  proper  form  in  several 

ways.  The  most  straightforward  procedure  is  to  bring  the  term 
2 2 

C / ( 1 + C ) d Tl  1 / dt  to  the  left-hand-side  and  simplify  to  obtain  the 
m m 1 r j 

following  expression: 


dt 
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(113) 


Another  way  of  handling  Equation  (112)  is  to  use  an  order  of  magnitude 

2 2 2 

argument  and  substitute  d T|. /dt  = - u;,  T),  on  the  right-hand-side.  Along 

. ^ ^ 10 
with  applying  the  MOA,  this  substitution  has  been  made  by  Culick  in 

deriving  Equations  (8.23)  and  (8.25).  With  this  substitution.  Equation 

(112)  becomes: 
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(114) 


It  is  now  apparent  that  the  right-hand-side  of  Equation  (114)  differs  from 
the  right-hand-side  of  Equation  (113)  by  the  factor  1/(1  + C ),  thus  it 
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is  expected  that  the  solutions  of  these  two  equations  will  differ  by  a 
similar  amount. 

Equations  (113)  and  (114)  respectively  describe  the  Galerkin  Equa- 
tions obtained  using  Culick's  formulation  without  and  with  the  above-men- 
tioned approximations.  The  MOA  can  now  be  applied  to  solve  each  of  these 
equations.  Letting 

Tlj^(t)  = g^(t)  sin  (ju^t)  + h^(t)  cos  (u)j^t) 

where  g^  and  hj^  are  slowly-varying  functions,  and  using  Equations  (32), 

equations  describing  the  behavior  of  gj^  and  hj^  can  be  derived.  Thus  with 

Culick's  formulation,  four  sets  of  results  were  obtained  for  a given  case: 

two  solutions  using  tne  Galerkin  method  both  with  (Equation  (114))  and 

2 2 2 

without  (Equation  (113))  the  substitution  d Tj  /dt  = ~ U)  ^ in 

and  the  two  corresponding  solutions  obtained  using  the  MOA. 

The  decay  rates  obtained  with  Culick's  formulation  for  cr  = 2 . 5 |j,  , 

C = 0.36,  and  f = 800  Hz  are  presented  in  Table  12  along  with  the  decay 
m e 

rates  obtained  with  the  present  analysis  (i.e.,  using  Equations  (25)  for 
the  Galerkin  solutions  and  Equations  (39)  for  the  MOA).  In  these  compu- 
tations the  thermal  energy  transfer  between  the  particles  and  the  gas  phase 
was  ignored.  For  this  case  the  "exact"  calculations  using  the  Kooker-Zinn 
model  yield  a = -57  sec 

A comparison  of  the  results  presented  in  Table  12  leads  to  several 

observations.  In  every  case,  the  MOA  yields  considerably  higher  decay 

2 2 2 

rates  than  the  Galerkin  method.  The  substitution  d T]j/dt  = - (t  results 
in  a significant  decrease  in  the  computed  decay  rates.  The  solution  pre- 
sented in  Reference  10  (i.e.,  with  MOA  and  making  the  substitution 
2 2 2 

d T|j^/dt  = agrees  well  with  numerical  results,  but  the  corres- 

ponding solution  using  the  Galerkin  method  does  not  agree  well.  Finally, 

the  Galerkin  equations  as  derived  by  the  present  analysis  and  that  of 

2 2 2 

Reference  10  without  the  substitution  d T]  j^/dt  = " uj  give  good  re- 

sults, but  in  both  cases  the  MOA  leads  to  incorrect  results. 

The  results  reported  above  raise  serious  questions  regarding  the 
applicability  of  the  MOA  over  the  wide  range  of  solid  rocket  operating 
conditions  experienced  in  practice.  The  example  considered  above  demon- 
strated that  good  results  could  be  obtained  with  the  MOA  in  a relatively 
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2 2 2 

simple  situation  only  after  the  approximation  d / dt  = " lU  j, was  intro- 
duced into  the  Galerkin  equation  while  erroneous  results  were  obtained 
when  the  MOA  was  applied  to  the  "exact"  Galerkin  equation.  At  this  point, 
in  view  of  the  absence  of  contrary  evidence,  one  cannot  help  the  feeling 
that  at  least  in  this  example  the  success  of  the  MOA  was  fortuitous. 


Table  12.  Comparison  of  Approximate  Solutions  for  Particles  in  a Box 


Galerkin 

Method  of 

Method 

Averaging 

Present  Approximate 
Analysis 

-55.7  sec"^ 

70  1 

-72,1  sec 

Analysis  of  Reference  (10) 
Without  the  Approximation 

2 

d ",  _ 2 ^ 

1 - -^1  'll 

dt 

-53.2  sec"^ 

-1  1 

-71.9  sec  1 

Analysis  of  Reference  (10) 
With  the  Approximation 

; ^ 1 'i 



-42,0  stc 

-52.8  sec  ^ 

The  applicability  of  the  MOA  in  the  analysis  of  solid  rocket  motors 
cannot  be  judged  solely  on  the  basis  of  the  above  comparisons  for  particles 
in  a box.  Therefore  more  extensive  comparisons  were  made  for  a hypotheti- 
cal motor  in  which  the  propellant  is  insensitive  to  the  pressure  oscil- 
lations (i.e.,  W = 0).  These  calculations,  therefore,  include  the  effects 
of  mean  flow,  flow  turning,  and  nozzle  damping  as  well  as  the  effects  of 
particle  damping.  Four  particle  diameters  were  considered,  2 p.  , 6 p.  , 10  p , 
and  20 p for  the  same  motor  geometry  and  steady-state  properties  consi- 
dered in  Section  4.2  (Figure  24).  Approximate  solutions  were  calculated 
using  the  following  three  techniques:  (1)  the  Galerkin  method  applied 

to  Equations  (13)  and  (14)  (i.e.,  numerical  solution  of  Equations  (22) 
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and  (23)),  (2)  the  Galerkin  method  applied  to  Equation  (17)  in  which  gas 
and  particle  equations  have  been  combined  (i.e.,  numerical  solution  of 
Equations  (25)),  and  (3)  the  MOA  applied  to  Equations  (25)  (i.e.,  numerical 
solution  of  Equations  (39)).  Decay  rates  and  frequencies  obtained  with 
the  three  approximate  techniques  are  compared  with  the  Kooker-Zinn  "exact" 
solutions  in  Table  13  for  3%  initial  disturbances. 

Several  interesting  observations  can  be  made  from  the  data  pre- 
sented in  Table  13.  For  the  larger  particle  diameters  of  lOp,  and  20p, 

(i.e.,  smaller  values  of  K) , all  three  approximate  solutions  are  in  good 
agreement  with  the  "exact"  solutions.  On  the  other  hand,  the  approximate 
solutions  for  2p,  and  bp,  diameter  particles  (large  K)  differ  from  the 
"exact"  solutions  by  various  amounts.  In  both  cases  the  Galerkin  methods 
yield  decay  rates  which  are  smaller  than  the  "exact"  decay  rates,  with 
Equations  (25)  yielding  the  greatest  error  (it  was  necessary  to  neglect 
the  spatial  derivatives  of  the  particle  velocity  in  order  to  combine  the 
gas  and  particle  equations).  The  MOA,  surprisingly,  yielded  decay  rates 
which  were  in  better  agreement  with  the  "exact"  decay  rates  than  the 
Galerkin  solutions,  in  spite  of  the  fact  that  additional  approximations 
were  made  in  applying  the  MOA  to  Equations  (25).  Apparently  the  error 
associated  with  applying  the  MOA  to  Equation  (25)  partially  compensated 
for  the  errors  associated  with  the  order  of  magnitude  approximations  and 
the  Galerkin  method  used  in  deriving  Equations  (25).  Since  it  is  impos- 
sible to  guarantee  that  this  fortuitous  circumstance  will  occur  in  every 
situation,  the  applicability  of  the  MOA  to  more  general  solid  rocket  in- 
stability problems  is  still  open  to  question. 

To  further  investigate  this  question,  a more  realistic  case  was 
considered  in  which  both  combustion  driving  and  particle  damping  were 
present.  Thus,  the  three  approximate  solutions  were  compared  with  the 
"exact"  solutions  for  the  case  shown  in  Figures  26  and  27  of  Section  4.2 
(i.e.,  a = 2.5p,,  = 0,36,  = 0.078,  = 4.17).  The  results  of  these 

calculations  are  shown  in  Figures  58  and  59.  Figure  58  shows  growth  rate 
versus  number  of  wave  cycles  for  3%  amplitude  IL  mode  initial  disturbances, 
while  Figure  59  shows  the  growth  toward  limiting  amplitude  for  10%  initial 
disturbances.  From  Figure  58  it  is  evident  that  the  Galerkin  method 
(Equations  (22)  and  (23))  underestimates  the  growth  rate  and  the  MOA  over- 
estimates the  growth  rate  by  about  the  same  amount  with  respect  to  the 
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Table  13.  Comparison  of  Approximate  Solutions  for  Motor  with  R = 0 


Particle  Size,  jx 

2.0 

6.0 

10.0 

20.0 

Drag  Constant,  K 

46.72 

5.191 

1.869 

0.4672 

Damping,  sec  ^ 

Gdierkin,  Eqs.  22 

and 

23 

-282 

-502 

-582 

-408 

Galerkin,  Eqs.  25 

-260 

-488 

-582 

-407 

M.O.A. , Eqs.  39 

-312 

-538 

-579 

-408 

Exact 

-322 

-525 

-589 

-411 

Frequency,  Hz 

Galerkin,  Eqs.  22 

and 

23 

975 

989 

1031 

1071 

Galerkin,  Eqs.  25 

975 

98  7 

1029 

1071 

M. O.A. , Eqs.  39 

979 

1000 

1041 

1074 

Exact 

976 

990 

1038 

1068 
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"exact"  solution,  while  the  Galerkin  method  (Equations  (25))  yields  growth 
rates  between  the  "exact"  and  MOA  values.  Figure  59  shows  that  the  MOA 
also  predicts  a considerably  higher  limiting  amplitude  than  either  the 
Galerkin  methods  or  the  "exact"  analysis.  These  results  show  that  the 
fortuitous  agreement  between  the  decay  rates  predicted  by  the  MOA  and  the 
"exact"  analysis  for  motors  without  combustion  driving  does  not  occur  when 
combustion  driving  is  present. 

In  conclusion,  the  considerable  reduction  in  computation  time 
obtained  with  the  MOA  clearly  warrants  its  use  in  the  nonlinear  analysis 
of  solid  rocket  motors  without  particles.  For  cases  in  which  particles 
are  present,  however,  the  MOA  should  be  used  with  caution,  particularly 
for  small  particles,  since  significant  errors  may  result. 

4,8  Computation  Time 

Table  14  shows  typical  computation  times  obtained  with  the  approxi- 
mate and  "exact"  programs  under  various  conditions  using  the  CDC  CYBER  70/74 
computer.  If  two  values  are  given  for  the  computation  time,  the  first  value 

refers  to  the  case  when  particles  are  absent.  These  results  show  that  for 

equal  time  step  sizes  the  two  Galerkin  approaches  are  about  equally  fast  and 
that  the  MCA  is  much  slower.  The  MOA,  however,  can  be  used  with  much  larger 

time  step  sizes  (one  or  two  steps  per  wave  cycle)  resulting  in  much  shorter 

computation  times  than  the  Galerkin  method.  The  T-burner  program  runs  about 
three  times  faster  than  the  corresponding  motor  program  due  to  the  simpler 
eigenfunctions  used  in  the  T-burner  analysis.  For  equal  time  step  sizes,  the 
approximate  programs  (excluding  the  MOA)  yield  from  a two-fold  to  a ten-fold 
increase  in  computational  speed  over  the  "exact"  programs.  For  cases  where 
the  MCA  yields  accurate  results,  computational  speeds  obtained  with  the  MOA 
may  approach  a hundred  times  that  obtained  with  the  "exact"  analysis. 


Table  14.  Comparison  of  Computation  Times 


Method 

Number  of 
Modes 

Time  Step 
Size 

Computation 

Time 

Seconds/wave  cycle 

Galerkin 

Eqs.  (22)  and  (23) 

4 

0.025 

7 

Galerkin 

Eqs.  (22)  and  (23) 

5 

0.025 

9-12 

Galerkin 

Eqs.  (22)  and  (23) 

6 

0.025 

15-21 

Galerkin  Eqs.  (25) 

6 

0.025 

20 

6 

0.025 

40 

MQA 

6 

0.25 

4 

T-Burner 

5 

0.025 

4 

Exact 

- 

0.06 

20 

Exact 

- 

0.025 

50 
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CONCLUSIONS  AND  RECOMMENDATIONS 


5,1  Conclusions 

The  primary  objective  of  this  research  program  was  to  evaluate  the 
usefulness  of  approximate  nonlinear  analysis  techniques  for  predicting  the 
stability  characteristics  of  solid  rocket  motors  and  T-burners.  This  work 
involved  the  adaptation  of  a previous  liquid-rocket  approximate  analysis 
for  the  solution  of  solid-rocket  instability  problems.  This  technique,  which 
was  based  on  the  solution  of  an  approximate  wave  equation  by  means  of  the 
Galerkin  method,  was  also  further  simplified  by  application  of  the  Method  of 
Averaging  (MQA),  The  usefulness  of  the  approximate  techniques  (Galerkin  and 
MQA)  in  the  nonlinear  stability  analysis  of  solid  rocket  motors  was  then  in- 
vestigated for  cases  in  which  gasdynamic  mode-coupling  was  the  only  nonlin- 
earity considered.  For  this  case  an  extensive  parametric  study  was  performed 
in  which  the  approximate  solutions  were  compared  with  available  "exact"  solu- 
tions, Later  the  effects  of  nonlinearities  in  the  combustion  driving  and  par- 
ticle damping  mechanisms  were  assessed  using  both  "exact"  and  approximate 
models,  which  required  the  development  of  heuristic  nonlinear  combustion  re- 
sponse and  particle  drag  laws  for  the  approximate  analysis.  Finally,  the 
applicability  of  the  approximate  model  to  the  nonlinear  analysis  of  T-burners 
was  determined. 

The  objective  of  this  program  has  been  reached.  The  predictions  of  the 
approximate  techniques  have  been  compared  with  the  "exact"  solutions  for  a wide 
variety  of  motor  operating  conditions.  On  the  basis  of  the  predicted  growth/decay 
rates,  frequencies,  limiting  amplitudes,  and  pressure  waveforms,  the  Galerkin 
technique  (without  averaging)  was  found  to  be  sufficiently  accurate  for  most  of 
the  cases  considered.  For  cases  in  which  particle  effects  were  important,  these 
studies  also  showed  that  the  MQA  was  generally  less  reliable  than  the  Galerkin 
method.  These  findings  were  also  supported  by  comparisons  with  experimental 
motor  data.  In  the  T-burner  studies,  it  was  found  that  the  approximate  T-burner 
analysis  failed  to  accurately  predict  limiting  amplitudes,  although  growth  rates 
at  low  amplitudes  were  in  fairly  good  agreement  with  measured  data. 

Additional  conclusions  reached  during  this  investigation  are  summarized 
below  in  four  parts:  accuracy  of  the  approximate  methods,  number  of  modes  re- 
quired to  satisfactorily  model  the  instability,  effects  of  the  improvements 
(i.e.,  incorporation  of  nonlinear  burning  rate  and  nonlinear  particle  drag  laws 
into  the  models),  and  computation  time. 
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Accuracy  of  the  Approximate  Method.  The  accuracy  of  the  Galerkin  method 
(without  averaging)  was  found  to  be  dependent  on  the  case  under  consideration. 
The  Galerkin  method  was  in  very  good  agreement  with  the  "exact"  solutions  for 
the  two  simplest  cases  considered;  (1)  a motor  without  particle  damping  or 
combustion  driving  (mean  flow,  flow  turning,  and  nozzle  damping  only)  and  (2) 
a particle-gas  mixture  in  a closed-ended  box.  For  a motor  with  particle  damp- 
ing but  without  combustion  driving,  the  approximate  analysis  underestimated 
the  damping  by  as  much  as  137o  for  particles  smaller  than  15  microns.  For  mo- 
tors with  combustion  driving  the  approximate  analysis  tended  to  overestimate 
both  the  growth  rates  at  low  amplitudes  and  the  limiting  amplitudes  for  cases 
without  particles,  while  these  quantities  were  generally  underestimated  for 
cases  with  particles.  For  these  cases  the  approximate  analysis  accurately 
predicted  the  expected  trends  regarding  the  effect  of  particle  size  and  con- 
centration and  propellant  lesponse  upon  limiting  amplitude.  For  T-burners  the 
approximate  analysis  also  coi  ectly  predicted  the  effect  of  various  particle 
properties,  propellant  response,  and  cup-grain  length,  but  for  given  values 
of  these  parameters  the  growth  rates  and  limiting  amplitudes  did  not  agree 
well  with  the  available  "exact"  solutions. 

The  accuracy  of  the  MCA  was  also  found  to  be  case  dependent.  For  small 
to  moderate  amplitude  disturbances  in  motors  without  particles  the  MOA  and 
Galerkin  method  were  equally  accurate.  For  small  particles  in  a box,  however, 
the  MCA  predicted  excessive  damping  compared  to  the  "exact"  analysis  and  the 
Galerkin  method.  Although  the  MCA  gave  better  agreement  with  the  "exact"  analy- 
sis than  the  Galerkin  method  for  a motor  with  small  particles  but  no  combustion 
driving,  the  MOA  did  not  yield  accurate  results  for  growth  rates  and  limiting 
amplitudes  in  the  same  motor  when  combustion  driving  was  Included, 

Comparison  of  the  approximate  (Galerkin)  solutions  with  available  motor 
test  data  gave  reasonable  agreement  with  the  measured  limiting  amplitudes  con- 
sidering the  difficulties  encountered  in  estimating  the  transient  combustion 
parameters.  The  approximate  analysis  correctly  predicted  the  experimentally 
observed  trend  of  increasing  limiting  amplitude  with  increasing  steady-state 
Mach  number.  For  T-burners  the  predicted  small  amplitude  growth/decay  rates 
were  in  good  agreement  with  the  measured  values  for  reasonable  values  of  the 
transient  burning  rate  parameters;  however,  the  predicted  limiting  amplitudes 
were  far  higher  than  the  measured  values. 
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Number  of  Modes  Required.  The  number  of  modes  required  to  adequately 
model  the  instability  depended  principally  upon  the  relative  stability  of 
the  higher  frequency  modes  with  respect  to  the  fundamental  mode,  A three  or 
four  mode  series  appeared  to  be  adequate  for  cases  in  which  the  higher  modes 
were  heavily  damped,  such  as  cases  in  which  small  particles  are  present.  When 
the  higher  modes  were  less  strongly  damped,  six  or  more  modes  were  necessary 
(i,e,,  motors  with  large  particles  or  without  particles).  The  most  unfavorable 
situations  arose  when  the  propellant  response  for  the  higher  modes  was  comparable 
to  or  greater  than  the  propellant  response  for  the  fundamental  mode;  here  a 
large  number  of  modes  were  required  and  the  approximate  solutions  required 
more  computer  time  than  the  "exact"  solutions. 

Effects  of  the  Improvements,  Results  obtained  with  the  "exact"  and  ap- 
proximate models  Indicate  that  combustion  nonlinearities  should  be  Included 
in  the  analysis  of  solid  rocket  instabilities  whenever  the  oscillation  ampli- 
tudes exceed  7-107,  and  sometimes  for  amplitudes  as  small  as  3%,  Depending  upon 
various  motor  and  propellant  parameters,  particularly  frequency,  combustion 
nonlinearities  may  either  decrease  or  increase  combustion  driving  with  increas- 
ing amplitude;  in  the  latter  case  pulsed  instabilities  may  result.  Although  the 
heuristic  nonlinear  combustion  model  developed  in  this  study  accounts  for  the 
above  effects,  it  needs  to  bo  improved  to  account  for  the  nonlinear  dependence 
of  the  combustion  driving  upon  amplitude  obtained  with  the  "exact"  analysis. 

On  the  basis  of  the  nonlinear  particle  damping  studies,  the  nonlinear 
particle  drag  law  should  be  used  for  the  higher  frequencies  and  larger  parti- 
cle sizes.  Depending  on  the  particle  size,  nonlinear  particle  drag  effects  may 
either  increase  or  decrease  particle  damping  with  increasing  amplitude;  in  the 
latter  case  pulsed  instabilities  may  result.  The  nonlinear  drag  effect  is  im- 
portant even  for  small  amplitude  disturbances  due  to  the  steady-state  contri- 
bution to  the  relative  veloci'  between  particles  and  gas.  Although  the  heuris- 
tic nonlinear  particle  drag  model  used  in  the  approximate  analysis  predicts  the 
above  trends,  more  parametric  studies  are  needed  in  order  to  guide  further  im- 
provements of  this  model. 


Computation  Time.  For  equal  integration  time-step  sizes,  the  Galerkin 
method  program  with  a six-mode  series  runs  about  twice  as  fast  as  the  "exact" 
program.  The  six-mode  series  requires  about  20  seconds  per  wave  cycle  central 
processor  time  on  the  CDC  CYBER  70/74  computer  used  in  these  studies.  For  cases 
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in  which  a four  mode  series  is  adequate,  a further  four-foid  increase  in  speed 
can  be  obtained.  The  MCA  program  runs  about  ten  times  faster  than  the  Galerkin 
program  because  a much  larger  time  increment  can  be  used.  Due  to  the  simpler 
eigenfunctions  used  to  represent  the  T-burner  solutions,  the  T-burner  program 
runs  about  three  times  faster  than  the  corresponding  motor  program. 

Under  the  most  favorable  circumstances  the  Galerkin  method  may  yield  a 
ten-fold  saving  in  computer  time  as  compared  to  the  "exact"  analysis,  while 
the  MCA  (when  applicable)  may  yield  a hundred-fold  increase  in  computational 
speed  over  the  "exact"  analysis.  Thus  the  approximate  solution  techniques  may 
be  used  for  extensive  parametric  studies  of  nonlinear  solid  rocket  and  T-burner 
stability  behavior  which  would  be  impractical  with  the  "exact"  analytical  tech- 
niques. 


5,2  Recommendations 

A number  of  investigations  aimed  at  improving  and  generalizing  the  ap- 
proximate models  are  recommended  for  future  research.  Some  of  these  can  be 
readily  accomplished  with  only  a modest  expenditure  of  time  and  effort,  while 
others  will  require  a major  research  effort  for  their  completion. 

One  of  the  basic  assumptions  used  in  the  development  of  the  present 
approximate  and  "exact"  analyses  was  that  the  exchange  of  thermal  energy  be- 
tween the  particles  and  gas  is  negligible.  Since  it  is  known  that  particle 
thermal  effects  may  contribute  significantly  to  the  total  particle  damping, 
they  should  be  incorporated  into  the  approximate  and  "exact"  models. 

The  present  approximate  program  for  motors  is  restricted  to  full-length, 
cylindrical  propellant  grains.  The  approximate  analysis  should  be  extended  to 
more  general  motor  geometries  such  as  partial  length  grains,  gaps  in  grain, 
end-burning  grains,  variable  port  area,  and  arbitrary  grain  cross-sections. 

Two  of  these  features,  partial  length  grains  and  end-burning  grains,  have 
already  been  incorporated  into  the  T-burner  model  and  could  easily  be  in- 
cluded in  the  analysis  of  motors. 

More  extensive  parametric  studies  are  needed  to  adequately  assess  the 
effects  of  nonlinear  combustion  driving  and  nonlinear  particle  damping  upon 
the  stability  characteristics  of  motors  and  T-burners,  These  studies  should 
be  conducted  using  both  the  "exact"  model  and  the  approximate  model.  The  re- 
sults of  these  studies  should  then  be  used  as  a guide  toward  improvements  of 
the  heuristic  nonlinear  combustion  response  and  particle  drag  laws. 
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In  view  of  the  limited  number  of  comparisons  with  experimental  data 
obtained  during  this  research  effort,  more  extensive  comparisons  between  the 
approximate  solutions  and  motor  test  data  are  badly  needed.  This  will  only 
be  possible  as  more  motor  test  data  becomes  available. 

More  studies  are  needed  to  determine  the  cause  of  the  poor  agreement 
between  T-bumer  limiting  amplitudes  predicted  by  the  approximate  model  and 
the  values  obtained  by  the  "exact"  analysis  or  from  T-burner  test  data.  This 
will  probably  involve  the  development  of  a nonlinear  vent  model  and  the  in- 
clusion of  wall  heat  transfer  effects. 

It  is  well  known  that  the  combustion  response  of  solid  propellants 
depends  upon  both  the  pressure  and  velocity  oscillations.  While  considerable 
work  has  been  done  to  date  on  the  investigation  of  the  pressure-coupled  re- 
sponse, relatively  little  work  has  been  done  on  the  elucidation  of  the  vel- 
ocity-coupled response.  It  is  therefore  recommended  that  the  available  in- 
formation on  the  velocity-coupled  response  be  used  to  develop  analytical 
expressions  capable  of  describing  the  velocity-coupled  response  of  solid 
propellants,  and  that  these  expressions  should  be  incorporated  into  the 
approximate  and  "exact"  nonlinear  instability  models.  Comparisons  should 
then  be  made  between  the  predictions  of  the  available  models,  with  and  with- 
out the  velocity-coupled  response.  Such  comparisons  should  be  useful  in  de- 
termining the  role  of  velocity-coupled  response  in  the  behavior  of  unstable 
rockets,  particularly  regarding  the  development  of  pulsed  instabilities. 

With  the  completion  of  the  above  recommended  improvements,  the  value 
of  the  approximate  analysis  as  a useful  engineering  tool  will  be  greatly 
enhanced. 
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APPENDIX  A 

DERIVATION  OF  APPROXIMATE  EQUATIONS 


A-1.  Derivation  of  Equations  (13)  and  (14). 

Under  the  assumption  stated  in  Section  2,  the  system  of  nondimens ional 
conservation  equations  that  describe  the  unsteady  behavior  of  the  combustor 
two-phase  flow  are  given  by  Equations  (1)  through  (6).  The  unsteady  flow  in 
the  T-burner  vent  region  can  also  be  described  by  these  equations  if  the  axial 
velocities  in  the  source  terms  arising  from  the  mean-flow/acoustic  interactions 
are  multiplied  by  the  factor  1 - (Section  2.6).  Thus  the  gas  and  particle 
momentum  equations  become; 
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Writing  the  energy  equation  in  terms  of  the  enthalpy,  Equation  (5)  yields: 
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Using  order  of  magnitude  approximations.  Equations  (1),  (2),  (A-1),  (A-2),  and 
(A-3)  will  now  be  combined  to  obtain  approximate  equations  governing  the  gas 
and  particle  velocity  potentials. 

Each  of  the  dependent  variables  is  first  expressed  as  the  sum  of  a 
steady-state  quantity  and  a perturbation  quantity,  thus 

P = 0 + p' 
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where  ( ) denotes  a steady-state  quantity  and  ( ) ' di  . ites  a perturbation 
quantity.  Equations  (A-4)  are  then  substituted  into  the  conservation  equations 
and  the  equation  of  state.  Subtracting  out  the  steady-state  equations  (obtained 
by  setting  the  time  derivatives  equal  to  zero  in  the  conservation  equations) 
yields  a system  of  equations  describing  the  perturbation  quantities.  These  equa- 
tions are  too  complicated  to  be  solved  by  the  Galerkln  method  and  must  be  further 
simplified  by  order  of  magnitude  approximations. 

It  is  now  assumed  that  each  perturbation  quantity  and  the  steady-state 
Mach  number  are  of  0(b)  where  • is  a small  ordering  parameter  that  is  a measure 
of  the  wave  amplitude.  Thus  the  perturbation  equations  are  simplified  by  neglect- 
ing all  terms  of  order  higher  than  e^,  which  includes  products  of  three  perturba- 
tions, products  of  the  mean  flow  Mach  number  (i.e,,  u)  and  two  perturbations, 

-2  -3 

products  of  u and  a perturbation  quantity,  and  terms  proportional  to  u . As  a 

result  of  the  small  Mach  number  assumption,  the  steady-state  pressure,  density, 
and  enthalpy  are  nearly  constant  throughout  the  chamber,  thus  p = 1,  p = 1,  and 
li  = 1 (the  error  Introduced  by  this  approximation  is  0(e:  )),  Furthermore  the 
steady-state  continuity  equations  for  the  gas  and  particles  show  that  the  steady- 
state  source  terms  m^  and  are  also  of  0(e),  while  the  linear  response  function 
(i.e,.  Equation  (9))  indicates  that  for  R of  0(1)  the  unsteady  source  terms  m^  and 

o ^ 

m^  are  of  O(s^),  Introducing  these  approximations  into  Equations  (1),  (2),  A(-l), 

P 

(A-2),  and  (A-3)  give: 
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while  the  equation  of  state  becomes; 


p'=  h'+p'+p'h' 


(A-10) 


In  order  to  further  simplify  these  equations,  the  following  velocity 
potentials  are  defined; 
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Further  order  of  magnitude  considerations  are  also  used  to  simplify  the  particle 
drag  terms  in  Equations  (A-7)  and  (A-8).  Here  it  is  assumed  that  the  particle 
loading  of  the  propellant  is  sufficiently  low  that  p^  is  0(e)  which  implies  that 
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is  0(e  ).  Since  the  relative  velocities  u-u  and  u'-u'  are  of  0(e)  the  second 
P P P 

and  third  particle  drag  terms  on  the  right-hand-sides  of  Equations  (A-7)  and 

(A-8)  are  of  0(e3)  and  can  be  neglected.  Introducing  these  approximations  and 
the  velocity  potentials  into  Equations  (A-5)  through  (A-9)  and  using  subscript 
notation  for  the  partial  derivatives  gives  the  following  second  order  conserva- 
tion equations: 
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In  order  to  combine  Equations  (A-12)  through  (A-16)  to  obtain  Equations 

13  14 

(13)  and  (14),  the  substitution  principle  used  by  Lighthill  and  Blackstock^ 
is  invoked.  This  principle  states  that  any  factor  of  a second  order  term  may  be 
replaced  by  its  equivalent  first  order  expression,  because  any  more  precise 
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substitution  would  result  in  the  appearance  of  higher  order  terms.  The  substi- 
tution principle  is  used  here  to  replace  p'  and  h^  in  the  second  order  terms 
with  equivalent  first  order  expressions  involving  $,  The  first  order  expressions 
which  are  needed  in  this  analysis  are  as  follows; 


p'  = 

h'  = -(y-1)$^ 
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which  are  derived  in  Reference  (23),  In  addition  the  steady-state  mass  source 

2m  /R  may  be  replaced  by  the  equivalent  quantity  du/dx  as  shown  by  the  steady- 
8 

state  continuity  equation.  Making  these  substitutions  into  Equations  (A-12), 
(A-14),  (A-16),  and  (A-10)  and  rearranging  terms  yields  the  following  expressions: 
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Equations  (A-18)  through  (A-21)  will  now  be  combined  to  obtain  a non- 
linear wave  equation  for  the  gas  velocity  potential.  Substituting  Equation 
(A-21)  into  (A-20)  and  collecting  terms,  the  energy  equation  becomes: 
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The  gas  momentum  equation  (i.e,,  Equation  (A-19))  is  integrated  with  respect  to 
X (the  constant  of  Integration  must  be  zero  according  to  the  argument  given  in 
Reference  (22))  to  give: 
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which  is  also  Equation  (63)  of  Section  2.6.  Setting  Q gives  the  corres- 

ponding equation  for  motors;  that  is.  Equation  (24)  of  Section  2,3.  To  proceed 
with  the  derivation  of  Equation  (13),  Equation  (A-23)  is  differentiated  with  re- 
spect to  t and  the  resulting  expression  is  substituted  for  (l/>)p'  in  Equation 
(A-22),  Furthermore  the  continuity  equation  (i,e.  Equation  (A-18))  is  substi- 
tuted for  p ' in  Equation  (A-22),  Combining  terms,  neglecting  the  nonlinear  par- 
ticle term  in  Equation  (A-22),  and  using  the  first  order  relation,  4 = 4 . 

tt  XX* 

yields 


tt 


2“^xt  + 


(>  + 1 


V ) ^ 4 
V dx  *t 


+ 24  4 + (v-l)4  4 

x xt  ' t > 


2ifi' 


+ Kpp(4^  - 4p^  ) 


(Y-l)PpUp  (4p)^j.  + 


R ^ 


(A- 24) 

If  2A'/R  is  replaced  by  ip',  Equation  (A-24)  becomes  Equation  (61)  of  Section 
2,6  which  describes  the  unsteady  flow  in  the  T-burner  vent  region.  For  motors 
the  source  term  in  Equation  (A-24)  is  related  to  the  pressure  perturbation 
through  the  linear  response  function  (i,e,,  Equation  (9)).  Substituting  Equa- 
tion (9)  into  the  source  term  and  using  ~ du/dx  and  the  first  order  re- 
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lation  p'  = - yields: 


R \ ’'^‘cdx^’t  (A-25) 

Substituting  Equation  (A-25)  into  Equation  (A-24)  and  setting  = 0 yields 
Equation  (13)  of  Section  2.2. 

The  particle  potential  equation  (i.e..  Equation  (14))  is  readily  derived 
from  the  particle  momentum  equation  (i.e.,  Equation  (A-15)).  Using  the  steady- 
state  particle  continuity  equation  to  replace  2m^/R  with  p du^/dx  ia  Equation 
(A-15),  dividing  by  p^,  and  rearranging  terms  gives; 


= - 

c dx  t 


(A-26) 


where  du  /dx  is  assumed  to  be  constant  (the  resulting  error  is  of  0(e^)).  As  in 
P 

the  case  of  the  gas  momentum  equation  (i.e..  Equation  (A-19))  Equation  (A-26) 
can  be  integrated  with  respect  to  x (the  same  argument  holds  regarding  the  con- 
stant of  integration).  Performing  this  integration  and  rearranging  terms  gives 
the  desired  particle  potential  equation: 


^^p^  ^p^^p^  ^ 2^%^x  = ■ V df  S 


(A-27) 


which  is  the  same  as  Equation  (62)  of  Section  2.6.  Setting  = 0 gives  Equation 
(14)  of  Section  2.2. 

A-2.  Derivation  of  Equation  (17). 


The  basic  assumption  used  in  deriving  Equation  (17)  is  that  the  terms 
involving  the  spatial  derivatives  of  the  particle  velocity  (both  steady-state 
and  perturbation)  are  negligible  compared  to  the  remaining  terms  in  the  par- 
ticle momentum  equation.  Neglecting  these  terms  as  well  as  those  involving  p' 
(since  p^  is  0(c^))  in  Equation  (A-8)  yields; 


+ Ku  = Ku 
P 


(A-28) 


Using  the  integrating  factor  e and  the  initial  condition  u'(0)  = u'  , the 

P Pq 

solution  to  Equation  (A-28)  is 


/ -Kt  Kt'  / ,./ 

u = Ke  e u dt 

P J 


(A-29) 


Since  the  second  term  decays  rapidly  as  time  increases,  it  can  be  neglected  to 
give; 


u'  (t)  = Ke 

P 


■Kt 


/ Kt' 

I u e dt 


(A-30) 


Introducing  the  velocity  potentials  given  by  Equations  (A-11),  Equation  (A-30) 
becomes : 


dx 


= Ke 


-Kt 


(A-31) 

To  derive  Equation  (17)  the  analytical  solution  for  $ given  by  Equation 
(A-31)  is  substituted  into  the  second  order  gas  momentum  equation  (i.e..  Equa- 
tion (A-19))  with  V^=  0 to  obtain; 


[ y-'-"  *t  * ‘ + i*] 


= 0 


(A-32) 


while  interchanging  the  order  of  differentiation  and  integration  in  the  last 
term  gives : 


^ r p 

dx  L Y 


+ +K»- * 


2r  „-Kt  r'" 


- K Pp  e 


J $e’^’''dt'  j = 0 

(A-33) 


Integrating  with  respect  to  x and  rearranging  terms  gives: 


■'  ■ - '[‘t 


2-  -Kt 
- K Pp  e 


f ^ -‘""'dt'  ] 


(A-34) 

which  is  the  same  as  Equation  (28)  in  Section  2.3. 

Next,  Equation  (A-31)  is  substituted  into  the  energy  equation  (i.e.. 
Equation  (A-22))  and  the  nonlinear  particle  term  is  neglected  to  obtain: 
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— - p*  - (v-1)  # $ - (v-1)  $ = 

Y t ’ t tt  ''  dx  t 


2ifi' 

__s 

R 


- 1)  - (Y  - 1)  ^ [p’papKe'^'^  J dt'] 


(A-35) 

As  in  the  previous  derivation  for  Equation  (A-24),  the  momentum  equation  (i.e,. 
Equation  (A-34))  is  differentiated  vith  respect  to  time,  which  yields: 

Y Pt  = - {^tt  ^ ^^xt  Vxt  - ^^t 


v2-  * . w-3-  -Kt  Kt'  / \ 

-Kp§+Kp  e $e  dt  ^ 

p p i ^ 


(A-36) 


Substituting  Equation  (A-36)  for  (1/y)p^  and  Equation  (A-18)  for  p'  in  Equation 
(A-35),  combining  terms,  and  using  the  first  order  relation,  4 = 4 . gives 

’ 1 1 XX ’ ^ 

the  desired  equation: 


du 


^xx  ■ ^t  = 2^'^xt  ^*t  + 


2- 


+ KPp$^  - (Y-1)  KPpUp4^  - K p 


+ p e"'^*^  4e*^'^  dt' 

P i 

^ / 1 \ ,/2-  - -Kt  , Kt'^  y 
+ (y-1)  K p u c 4 e dt 

p p J X 


2ifi' 


R ^c 


(A-37) 


Finally  replacing  2i?iVr  with  Equation  (A-25)  yields  Equation  (17)  of  Section 

2.2. 
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APPENDIX  B 


USE  OF  COMPLEX  VARIABLES  IN  THE  SOLUTION  OF 
NONLINEAR  DIFFERENTIAL  EQUATIONS 

It  is  often  convenient  to  use  complex  variables  in  the  solution  of  the 
linear  equations  which  arise  in  acoustics,  combustion  instability  and  related 
fields.  In  this  case  the  solution  is  expressed  in  complex  form,  and  the  leal 
part  represents  the  physically  meaningful  solution.  However,  care  must  be  used 
when  applying  this  technique  in  the  solution  of  nonlinear  equations.  The  dif- 
ficulties that  are  encountered  in  applying  the  complex  variable  technique  to 
nonlinear  problems  will  be  illustrated  by  analyzing  the  following  simplified 
example.  Consider  the  nonlinear  wave  equation  given  by: 


$ 


XX 


\t 


1$ 


'■i-D 


A complex  solution  of  Equation  (B-1)  of  the  form  $ = cp  + iY  would  be  useful 
only  if  its  real  part,  cp,  satisfies  Equation  (B-1),  which  would  be  the  case  if 
the  equation  were  linear.  However,  straightforward  substitution  of  $ = cp  + iY 
into  Equation  (B-1)  and  separating  its  real  and  imaginary  parts  yields  the 
following  equation  for  cp : 


cp  - cp  = cpc?  - YY 
^xx  ^tt  ^t  t 


(B-2) 


indicating  that  the  real  part,  cp,  does  not  satisfy  Equation  (B-1)  because  of 
the  extra  term,  -YY^.,  appearing  on  the  right  hand  side.  In  order  to  eliminate 
this  extra  term,  the  form  of  the  original  differential  equation  (i,e,.  Equation 
(B-1)  must  be  modified. 

Since  Equation  (B-1)  supposedly  describes  some  physical  phenomenon,  and 
since  only  the  real  part  of  the  ccxnplex  solution  is  physically  meaningful,  then 
the  nonlinear  term  should  really  be  expressed  as  the  product  Re($)  x Re($  ) 
which  is  equivalent  to  (^ij^  + + $ $^  + $ $^  ) /4.  Substituting  this  expression 

into  Equation  (B-1)  yields; 

i - + $*$*1/4  (B-3) 

xxtt‘-t  t t t-' 
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Substituting  $ = cp  + iY  into  Equation  (B-3)  and  separating  its  real  and  imagi- 
nary parts  yield: 


cp  - cp 
XX  ^tt 


W, 


Y - Y =0 

XX  tt 


fB-4) 


which  shows  that  the  real  part  of  the  solution  of  Equation  (B-3)  satisfies  the 
desired  equation  (i.e..  Equation  (B-1))  and  the  imaginary  part  satisfies  a 
homogeneous  linear  wave  equation.  This  technique  was  applied  to  the  solution  of 
nonlinear  combustion  instability  problems  (i.e.,  to  Equation  (13)),  and  the  re- 
sulting modified  wave  equation  was  solved  using  the  Galerkin  Method,  Due  to  the 
approximate  nature  of  the  Galerkin  Method,  however,  the  resulting  solution  con- 
tained an  error  term  which  grew  without  limit.  Consequently,  the  above  procedure 
had  to  be  modified  in  order  to  obtain  satisfactory  solutions  of  Equation  (13) 
using  the  Galerkin  Method. 

An  alternate  technique  is  to  modify  Equation  (B-1)  such  that  both  the 
real  and  imaginary  parts  satisfy  the  original  equation.  This  can  be  done  by  re- 
placing terms  of  the  form  with  Re($)Re($^)  + ilm($)lm($^) ; using  the  rela- 


tions : 


Re(J)Re($ 


t>  ■ ( H-^)(  -h-^)  ■ »*<  ] 


/4 


(B-5) 


i 


ilm($)lm($^)  = -i  ~^~2 — " 


* * * * 


]« 


in  Equation  (B-1)  gives: 

^xx“  ^tt  ""  ^ ^ $j-)]  /4 


(B-6) 


Substituting  $ = cp  + iY  into  Equation  (B-6)  and  separating  into  its  real  and 
imaginary  parts  gives: 


Cp  - cp  = CPCP 
^xx  ^tt  ^t 


Y - Y = YY 
XX  tt  t 


(B-7) 
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which  shows  that  both  cp  and  Y satisfy  Equation  (B-1).  Applying  this  method  to 
the  solution  of  Equation  (13)  yields  the  following  modified  wave  equation: 


'xx-  ^t  = 2“^t  ^ g - vKh^  ^ 


[v  (Y-i)PpUp(y^, 


+ $%*  1 + i±i  ] 

2 L X xt  X xt  J 2 L X xt  x^xt  J 


{(1-i)  ] + (1+i)  [^C  + ^t^xx]} 


Application  of  the  Galerkin  Method  to  Equation  (B-8)  yields  Equations  (22)  of 
Section  2,3 
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APPENDIX  C 


COEFFICIENTS  APPEARING 
MODE  AMPLITUDE 


IN  THE  APPROXIMATE 
EQUATIONS 


C-1.  Coefficients  Appearing  in  Equations  (22)  and  (23). 

The  coefficients  of  the  linear  terms  in  Equations  (22)  are: 
1 

C (i,m)  = f X X* 
o^-"  J m j 


dx 


for  1 S m S 1 


C^(j,m)  = 


C^ ( j ,m) 


-J 


-1  d^X 


m * 

® X.  dx 


o dx 
0 


2 j 


* 


for  N+1  S m s 


for  1 S m S 1 


for  N+1  ^ m £ 


dX 


C2(j»™)  = 


"I  3^  ’‘j 


+ (v+1)  f ^ X X*  dx 
' ' J dx  m J 


+ yY  X^(1)  X.(1) 


for  I ^ m s 1 


^2(1."!)  = 


- 


dX 


-(Y-l)PpJ  Up(x)  — X.  dx 


for  N+1  Sms 


C3(j,m)  = p 


* 

X X. 

1 J m J 


dx 


for  1 5 m ^ r 


l(j.m)  = -p'pf  X^X*  d> 


for  N+1  S m < 
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(C-1) 

2N 

(C-2) 

2N 

(C-3) 

2N 

(C-4) 

2N 


I 


* 


n 


X- 
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for  1 S m S N 


S«'“>  ■ -'fsVj'"' 
= 0 


(C-5) 

for  N+1  S m S 2N 


The  coefficients  of  the  linear  terms  in  Equations  (23)  are: 


C^Cj.m)  = 
CgCj.m)  = 
Cg(j,m)  = 
C^(j,m)  = 


for  1 ^ m ^ N 


(C-6) 


* 

X X. 

J m j 

* 

-K  X X.  dx 

J ml 

o 

f u (x)  ^ X.  dx  + 1 X X.  dx 

J p dxj  Jdxmj 


for  N+1  m ^ 2N 


for  1 S m S N 


»1  du_ 


(C-7) 


+ K 


pi  * 
; X X. 
J m J 


dx 


for  N+1  S m S 2N 


The  coefficients  of  the  nonlinear  terms  in  Equations  (22)  are  given  below 
for  Isj^N,  l^m^N,  and  1 S n S N,  For  all  other  values  of  j,  m,  and  n, 

^ = °2  = °3  = ^ = °- 


1+i 

2 


1 dX 
, p m 

dX 

J dx 
o 

dx 

(Y-l)(l- 

4 

pi 

. 1 "> 

dx' 

] 

J dx 
o 

dx 

(v-D(i-n) 

dx 


pi  * 

f — ^ X X dx 

i dx^  " J 


(C-8) 


.1  d^X 


m 


dx 


* * 

X X,  dx 

n j 


(C-9) 
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D3(j ,m,n) 


1 -1  dX  dX 

l-l-i  p m n * 

2 J dx  dx  j 


dx 


2 * 

d X * 

m " 

X X,  dx 

n j 


1-i 

2 


J 

2 

o 

dx 

* 

n 

★ 

X.  dx 
J 

•y  * 
d^X 

m 

c dx 


X X,  dx 

n j 


(C-10) 


(C-11) 


In  the  integrals  appearing  in 
functions  X^ (x)  are  defined  by  Equat 
conjugates. 


the  coefficients  given  above,  the  eigen- 
★ 

ion  (20)  and  X^ (x)  denotes  their  complex 


C-2.  Coefficients  Appearing  in  Equations  (25). 

The  linear  coefficients  in  Equations  (25)  are  given  below  for  1 S j s N 
and  1 S m S N,  The  nonlinear  coefficients  in  Equations  (25)  are  identical  to 
those  in  Equations  (22)  and  are  given  by  Equations  (C-8)  through  (C-11). 


Cg(j,tn)  = 

X X.  dx 
J m J 

(C-12) 

fl  * 

2 if 

C^(j,m)  = 

- f — X . dx  - 

i dx2  J 

K.  p X X.  dx 

p Jq  ® j 

- k(y-Dp  j 0 (X) 
^ 0 ^ 

dX  ^ dX  ^ 

r-Sl  X.  dx  + ~ (l)X.(l) 

dx  J dx  J ' ' 

(C-13) 

C2(j,m)  = 

r^.  * 

2 u(x)  — X.  dx 

+ (y+1)  f ^ X X*  dx 

^ J dx  m J 

+ Kp  X X*  dx  + 

P ^ m J 

yYX  (1)X*(1) 

j 

(C-14) 
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C3(j,m)  = 


(C-15) 


X X*  dx 
tn  J 


+ 


, pi  dX 

Up(x)  dx 


c^a.m)  = 


X X*  dx 
m j 


(C-16) 


C-3,  MethodoloRv  for  Calculating  Coefficients  In  Equations  (39), 

Separation  of  Real  and  Imaginary  Parts.  Equations  (39)  were  derived  by 
applying  the  Method  of  Averaging  (MOA)  to  Equations  (25).  This  procedure  re- 
quires that  Equations  (25)  be  separated  into  their  real  and  imaginary  components. 
Assuming  that  A.(t)  = Fj(t)  + iGj(t),  substituting  into  Equations  (25),  and  sep- 
arating real  and  imaginary  parts  yields: 

dB 
m 

dt 


2N 

z 


C'(j,m) 


d^B 


dt 


C^'(j,m)B^+  C'(j,m)  +h^R^E'(j,m)  +h^R^E2(j,m) 


2N  2N 


+ C'(j,m)e"'^‘'  pB^(t')e’^*''dt'  + X!  S! 

o m=l  n=l 


dB 

D'(j.m,n)Bm^  - 0 


j = 1,2,3  ...  2N  (C-17) 

where  the  B^s  are  related  to  the  F^s  and  G^s  by  Equations  (33)  of  Section  2.4. 
m mm 

The  real  coefficients  C^,  C^,  C^,  C^,  E^,  and  d'  in  Equation  (C-17)  are  re- 
lated to  the  original  complex  coefficients  (i.e.,  C^,...C^,  D3,...D^)  appearing 
in  Equations  (25)  as  follows: 


Cj^(2j-1,  2m-l)  = 

Re[Cj^(j,m)] 

C^(2j-1,  2m) 

-Im[Cj^(j,m)] 

C^(2j,  2m-l) 

Im[Cj^(j,m)] 

C^(2j,  2m) 

Re[Cj^(j,m)J 
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for  k = j — N,  in  — N, 


E'(2j-1,  2m-l)  = Re[c^(j,m)] 


E'(2j-1,  2m)  = -Im[c^(j,m)J 


e'(2j,  2m-l)  = Im|c^(j,m)J 


E'(2j,  2m) 


[C4(j,m)] 


E'(2j-1,  2m-l)  = 


‘[C4(j,m)] 


E'(2j-1,  2m) 


•[C4(j,m)] 


e;(2J,  2m-l) 


•[C4(j,m)] 


^2(23,  2m  ) 


'[C4(j,ni)] 


(C-19) 


for  j = 1,2...  N,  m = 1,2,...  N,  and; 


D'(2j-1,  2m-l,  2n-l)  = Re[D^(j,m,n)  + D2(j,m,n)  + D3(j,m,n)  + D^(j,m,n)] 
D'(2j-1,  2m-l,  2n)  = Im[-D^  ( j ,m,n)  + D2(j,m,n)  - D^Cj.m.n)  + D^(j,m,n)J 

D'(2j-1,  2m,  2n-l)  = Im  ( j ,m,n)  - D2(j,m,n)  + D3(j,m,n)  + D^(j,m,n)J 

D'(2j-1,  2m,  2n)  = Re [-D^ ( j ,m,n)  + D2(j,m,n)  +D3(j,m,n)  - D^(j,m,n)] 


D'(2j,  2m-l,  2n-l)  = Im|p^(j,m,n)  +D2(j,m,n)  +D2(j,m,n)  +D^(j,m,n)j 

D'(2j,  2m-l,  2n)  « Re|p^(j,m,n)  - D2(j,m,n)  + D3(j,m,n)  - D^(j,m,n)] 


I 


d'(2j,  2m,  2n-l) 


d'(2j,  2m,  2n) 


Re[Di(j,m,n)  + D2(j,m,n)  - D2(j,m,n)  - D^(j,m,n)J 
Imj|-D^(j,m,n)  + D2(j,m,n)  + D2(j,m,n)  - D^(J,m,n)J 


(C-20) 


for  j = 1,2,...  N,  m = 1,2,...  N,  n = 1,2...  N. 


Transformation  to  Uncoupled  Form.  Equations  (C-17)  are  coupled  in  the 
second  derivatives;  that  is,  there  are  two  or  more  terms  in  each  equation. 
This  coupling  results  from  the  non-orthogonality  of  the  axial  eigenfunctions. 
In  order  to  apply  the  MOA  to  Equations  (C-17),  they  must  be  decoupled  to  the 
form : 


d^B, 


= f 


dt 


/ ^ ^ ^2n\ 

B2,...  B2j^,  I 


(C-21) 


in  which  only  one  second  derivative  appears  in  each  equation.  Using  Equation 
(C-21),  it  is  seen  that  Equation  (34)  (or  Equation  (C-17))  can  be  expressed  as 


c;f  = g 


(C-22) 


vdiere  C^  is  the  2N  x 2N  matrix  of  coefficients  of  the  coupled  system,  f is  the 


column  matrix  corresponding  to  the  right -hand -side  of  the  decoupled  system,  and 
g is  the  column  matrix  corresponding  to  the  right-hand-side  of  the  coupled  sys- 
tem. To  decouple  Equations  (C-17),  therefore.  Equation  (C-22)  is  solved  for  f, 
thus 

(C-23) 


where  C^^  is  the  inverse  of  the  matrix  C^.  Performing  these  operations  and 
equating  the  coefficients  of  like  terms  in  f and  C^^  g gives  the  following 


relations ; 


2N 


= Z!  c“^(j,k)  cf(k,m)  i = 1,2,3 


k=l 


“S’ (jin') 


'Et  (j,m,n)  = ^ D'(k,m,n) 

k=l 

where  "c^,  and  "d  are  the  corresponding  coefficients  of  the  decoupled  system. 

Due  to  the  complicated  nature  of  the  relationships  given  by  Equations 
(C-24),  explicit  expressions  for  the  coefficients  7^,  and  D in  terms  of  the 
spatial  integrals  can  not  be  given.  The  coefficients  must  be  calculated  numer- 
ically by  first  computing  the  complex  coefficients  using  Equations  (C-12)  through 
(C-16)  and  Equations  (C-8)  through  (C-11),  next  calculating  the  coefficients  of 
the  equivalent  coupled  real  system  using  Equations  (C-18),  (C-19),  and  (C-20), 
and  finally  computing  the  coefficients  of  the  decoupled  system  using  Equations 
(C-24).  A similar  procedure  is  used  in  the  numerical  solution  of  the  equations 
obtained  using  the  Galerkin  method  without  averaging  (i.e..  Equations  (22)  and 
(23)  or  Equations  (25)). 

Application  of  MQA.  The  MOA  procedure  is  described  in  Section  2.4  which 
results  in  Equations  (39).  In  performing  the  time  integrations  indicated  by 
Equations  (32),  two  integrals  of  the  products  of  three  trigonometric  functions 
arise.  These  appear  in  Equations  (39)  and  are  defined  by: 


2tt/ou, 

/%  J. 

<"1  J 

cos(u)  t) 

P 

COS 

o 

tt/2  if 

p = q 

+ Z- 

or 

q = p 

+ 't 

or 

-t  = p 

+ q 

0 for  all  other  integer  values  of  p,q,'t 
2n/u)^ 


= ui 


cos(u)  t)  sin(u)  t)  sin(uj,t)  dt 
1 J p q -t 


n/2  if  I = p + q 

or  q = p + -t 

-tt/2  if  p = q + <, 

0 for  all  other  integer  values  of  Pjq.-t 


(C-25) 


(C-26) 
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C-4,  Mode -Amplitude  Equations  for  T-Burners. 

The  mode -amplitude  equations  for  T-burners  are  similar  to  Equations  (22) 
and  (23)  with  an  additional  equation  for  u^  (Equation  (76)),  The  T-burner 
equations  are: 


2N  / d^A 

X)  1 ~~T  KC3(j,m) 

m=l  I L 


+ h^ft„C^(j.m)J  C^(j) 


2N  2N 


ZX  ''rndT  + ""m  dl 


m=l  n=l 


n , n 


dA  -k  ^ 

+ D2(j,m,n)  A^  ^ \ dT 


j = 1,2...  N 


(C-27) 


2N  / dA  \ 

m=l  ' J 


j = 1,2,. 


(C-28) 


uu  a OA 

dT  ^ <V^eff  ><  = <v/^eff  >X  df 


(C-29) 


The  linear  coefficients  appearing  in  Equations  (C-27)  are: 


Cjj(j.m) 


r X X 
J m j 


for  1 £ m & N 


(C-30) 
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for  N+1  S m £ 2N 


1 d^X 


. -J 


for  1 ^ m ^ N 


Cj (j,m) 


for  N+1  S m i 2N 


ri 

C2(j,m)  = 2j  G(x)  Xj  dx 


1 . (l-^^^)/2 

+ (y+1)  [ ^ X X dx  - V.  [ du/dx  X X.  dx 

Jaxmj  -oJ 


'=b[‘  + V'> 


(i-e^)/2 


for  1 S m 5 N 


-1  . _ dX 

C,(j,ni)  = -(\-l)  u p - — X.  dx  for 

2 J p dx  j 


N+1  S in  S 2N 


for  1 S m S N 


-f  p X X.  dx 
J p m j 


for  N+1  S m S 2N 


f X X.  dx  + f ^ X X.  dx 

J dx  m J J dx  m J 
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""b  [‘ 


for  1 s:  m s N 


for  Nfl  s m S 2N 


^ r- 
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(C-36) 


C^(j.n>) 


for  1 S m <;  N 


X X.  dx 
J m j 

0 

1 dX  dfl 

1 u (x)  3—^  X.  dx  + 3—^  X X.  dx 

J “p''  ''  dx  3 J dx  m J 


_ V r 3-^  X X.  dx  + K X X. 

V.J  dx  m J J m J 


(C-37) 


dx 


(l-P^)/2 


for  N+1  < m 5 2N 


The  nonlinear  coefficients  appearing  in  Equations  (C-27)  are  similar 
to  those  appearing  in  Equations  (22)  and  are  given  by. 


_ .1  dX  dX 

DjCj.m.n)  = J d^ 
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^ (Y-l)Clzil  15"  X X.  dx 

^ i dx2  ^ J 
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APPENDIX  D 


DERIVATION  OF  STEADY-STATE  SOLUTIONS 

D-1.  Steady-State  Solutions  for  Motors. 

The  steady-state  solutions  for  gas  and  particle  velocity  and  particle 

density  for  motors  with  full-length,  cy lindrtcally  perforated  grains  (i.e,, 

Equations  (40),  (41),  and  (42)  of  Section  2.5)  will  now  be  derived.  For  small 

steady-state  Mach  numbers  the  gas  density  and  mass  burning  rate  are  assumed 

to  be  constant  (the  variation  in  p and  m is  of  0(M^));  thus  for  p(x)  = 1 the 

g e 

steady-state  continuity  equation  becomes: 


du 

dx 


2ifi 

R 


(D-1) 


Integrating  Equation  (D-1)  and  using  the  head-end  boundary  condition  u(0)  = 0 


gives ; 


'2ili 


u(x)  = 


X = U X 

e 


(D-2) 


which  is  Equation  (40)  of  Section  2.5. 

The  steady-state  particle  density  and  velocity  are  obtained  from  the 

steady-state  particle  continuity  and  momentum  equations.  Neglecting  variations 

-2  - 

of  Pp  with  X as  being  of  O(M^),  dp^/Sx  5®  0 and  Equations  (2)  and  (4)  of  Section 
2.1  become: 


du 

E 

dx 


A 


(D-3) 


du 


P u j 
^ p p dx 


= p K(u  - u ) - ^ m u 
p-”  R p p 


(D-4) 


Integrating  Equation  (D-3)  with  Up(0)  = 0 gives: 


2m 

1 Hp  , 


(D-5) 


To  determine  p^,  Equations  (D-2),  (D-3)  and  (D-5)  are  substituted  into  Equation 
(D-4)  to  obtain  the  following  quadratic  equation  for  p : 


188 


* 


A 


2 


2 


0 


(D-6) 


-2 


P 


P 


which  is  independent  of  x, 
y ie Ids : 


P 


P 


Solving  Equation  (D-6)  and  taking  the  positive  root 


(D-7) 


which  is  Equation  (42)  of  Section  2.5.  The  particle  velocit”  is  then  obtained 
by  substituting  Equation  (D-7)  for  p in  Equation  (D-5)  to  obtain: 


Up(x)  = 


2u  X 
e 


I + All  + 


8Uf 


which  is  Equation  (41)  of  Section  2.5. 


(D-8) 


D-2.  Steady-State  Solutions  for  T-Burners. 

The  steady-state  solutions  for  gas  ana  particle  velocities  and  particle 
density  for  a cup-grain  T-burner  will  now  be  derived.  Solutions  are  obtained 
for  Regions  1,2,  and  3;  the  solutions  in  Regions  4 and  5 are  readily  obtained 
by  symmetry. 


Region  1 . The  steady-state  solutions  in  the  region  of  lateral  surface 
burning  (i.e.,  0 < x < 3/2  ) for  the  T-burner  differ  from  the  solutions  given 
previously  in  Section  D-1  for  motors  only  if  end-burning  is  present.  For  this 
case  the  steady-state  velocity  of  the  gases  leaving  both  the  lateral  propellant 
surfaces  and  the  end  surfaces  is  denoted  by  Uj^ ; hence,  m^  = u^^  for  p = 1.  In- 
tegrating the  steady-state  continuity  equation  (i.e..  Equation  (D-1))  and  using 
the  boundary  condition  u(0)  = u^  gives; 


u(x) 


u^(  1 + ^ X ) 


(0  X < 3/2)  (D-9) 


which  is 
As 

partic le 


Equation  (48)  of  Section  2.6. 

before,  the  particle  properties  are  obtained  by  solving  the  steady-state 
continuity  and  momentum  equations.  If  the  approach  used  in  Section  D-1  is 
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employed,  the  resulting  solution  for  is  dependent  on  x which  violates  the 
assumption  that  is  constant.  Furthermore  p -•  co  as  x -•  0 which  is  physi- 
cally meaningless.  Thus  for  the  case  of  end  burning,  the  assumption  of  a 
uniform  particle  density  is  inconsistent  with  the  particle  continuity  and 
momentum  equations. 

To  obtain  solutions  for  Pp(x)  and  Up(x)  an  alternative  solution  pro- 
cedure is  used.  Allowing  Pp  to  vary  but  assuming  that  is  a constant,  the 
particle  continuity  equation  becomes: 


( p u ) = 

P 


(D-10) 


Using  the  boundary  condition,  pp(0)  Up(0)  = , Equation  (D-10)  is  integrated 

to  obtain: 


= *p<  1 ^ R - > 


(D-11) 


The  steady-state  particle  velocity  is  then  obtained  from  Equation  (D-4)  by  sub- 
stituting Equation  (D-11)  for  PpOp  , substituting  Equation  (D-9)  for  u,  and 
approximating  Pp  in  the  particle  drag  term  by  = m^/  Simplifying  the  re- 

sulting expression  yields  a linear  first  order  differential  equation  of  the 


u = K. 
P 


(D-12) 


where  B = K/u,  + 2/R  and  A = 2/R.  This  equation  has  the  general  solution: 
0 


Gp(x)  = (1  + Ax  )'®^^  ( 1 + Ax  ) 


(D-13) 


Using  the  boundary  condition  Up(0)  = u^^  (the  particles  and  gas  are  assumed  to 
leave  the  propellant  surface  with  the  same  velocity),  introducing  the  expressions 
for  A and  B,  and  introducing  T]  = 4u^^/RK,  Equation  (D-13)  becomes: 


_ TJ+2 


(D-14) 


The  steady-state  solution  for  particle  density  is  obtained  by  substituting 


'•  * —v  j- 
. * T-  r 


Equation  (D-14)  into  Equation  (D-11)  and  solving  for  p^,  thus; 


Pp(x)  = 


1 + T1 


(D-15) 


1 -M  [1+  I X ] 


■(2T>4-2)/T1 


Equations  (D-14)  and  (D-15)  correspond  to  Equations  (49)  and  (50)  of  Section 
2.0  respectively 

Region  2.  In  Region  2 the  source  terms  m and  m vanish  and  the  steady- 

8 P 

state  continuity  equation  for  the  gas  phase  becranes; 


^ (pu  ) = 0 (D-16) 

which  is  integrated  to  obtain  pu  = constant.  The  constant  is  obtained  by  match- 
ing the  solutions  at  the  boundary  between  Regions  1 and  2.  Evaluating  Equation 
(D-9)  at  X = p/2  yields  the  constant  value  of  u in  Region  2,  thus: 

G(x)  = 0^(1  + p/R  ) (p/2  < X < (l-e^)/2)  (D-17) 


which  is  Equation  (55)  of  Section  2.6, 

The  particle  conservation  equations  in  Region  2 become: 


d 

dx 


P 

P P 


) = 0 


(D-18) 


- . 

*^p^p  dx 


p K(u-u  ) 
P P 


(D-19) 


Integrating  Equation  (D-18)  and  matching  solutions  at  the  boundary  between 
Regions  1 and  2 gives : 

p a = p'  (P/2)  G^(3/2)  (D-20) 

P P P P 

Substituting  Equation  (D-20)  into  Equation  (D-19),  approximating  p^  by  in 
the  right-hand-side  of  Equation  (D-19),  and  substituting  u = u(P/2)  yields  a 
linear  first  order  differential  equation  of  the  form: 
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where 


+ Au  = Q 
dx  p 


A = C K/[5  0/2)  u (e/2)] 
m ^ p p 


Q = KG(e/2)/[PpO/2)  Up(e/2)] 


Equation  (D-21)  has  the  general  solution  given  by: 


u (x)  = R/A  + C,  e 

P J- 


(D-21) 


(D-22) 


Matching  particle  velocities  at  the  boundary  between  Regions  1 and  2 gives  the 
constant  , Introducing  the  values  of  , A,  and  Q into  Equation  (D-22)  gives 
the  desired  steady-state  solution  for  the  particle  velocity  in  Region  2: 


u (x)  = u (P/2)  {l  - Y 


-o(x-3/2) 


(P/2  < X < (1  - p^)/2) 


(D-23) 


where  a = K/uO/2)  and  1]  is  the  same  as  in  Region  1,  Equation  (D-23)  corresponds 
to  Equation  (56)  of  Section  2.6. 

Substituting  Equation  (D-23)  for  u in  Equation  (D-20)  and  solving  for  p 


yields ; 


Pp(x)  = 


Pp(e/2)Up(p/2) 

G(e72) 


I 1 - 
1 1 + 


-o'(x-p/2) 


(D-24) 


(p/2  < X < (l-p^)/2) 


Region  3.  In  the  vent  region  the  source  terms  appearing  in  the  continuity 
and  momentum  equations  are  determined  by  the  flow  out  the  vent.  Assuming  uniform 
flow  across  the  vent,  the  steady-state  continuity  equation  for  the  gas  phase 


becomes ; 


(pu)  = - 2pG^(  1 + e/R)/e^ 


(D-25) 
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Again  assuming  that  p = 1 Equation  (U-25)  is  integrated  to  obtain 


u(x)  = G.(l  + e/R)(l  - 2x)/B 
b V 


(d-e^)/2  < X < (1  + e^)/2) 


(D-26) 


where  the  constant  of  integration  was  obtained  by  matching  the  solutions  at 
the  boundary  between  Regions  2 and  3,  Equation  (D-26)  is  the  same  as  Equation 
(70)  of  Section  2.6, 

Assuming  that  the  particle  sink  at  the  vent  is  simply  the  product  of  C 

m 

and  the  gas  sink,  the  particle  conservation  equations  become: 


— (p  G ) = - 2 C G,  ( 1 + p/R)/e 

dx  '^p  p'  mb' 


(D-27) 


P u = p K(u-u  ) + 2 C u,  ( 1+3/ R)  (1-V,)u  /p 

^p  p dx  ^p  ' p^  m b'  \ p/s-y  (D-28) 


Following  the  method  used  in  Region  1,  Equation  (D-27)  is  integrated  to  obtain: 


»p“p  • 


(D-29) 


which  is  substituted  into  Equation  (D-28)  to  obtain  a linear  first  order  dif- 
ferential equation  of  the  same  form  as  Equation  (D-12)  where  A = -2  and: 


® Gj^(l+8/R)  ■ 


Substituting  these  values  of  A and  B into  the  general  solution  given  by  Equation 
(D-13)  and  matching  solutions  at  x = (1-0  )/2  gives  the  desired  solution  for  u : 

V p 


Up(x)  = 


u^^(H-0/R) 


r i-2x  „ ri-2x  n'’  i 


(D-30) 


(4^ 


(1+0/R) 


(l-e^)/2  < X < 1/2 


• _v  r, 

A • ^ : 


and 


q = 2/n^  + (V_^-l) 


Substituting  Equation  (D-30)  Into  Equation  (D-29)  and  solving  for  gives: 


Pp(x)  = 


‘-u^rT 


(1  - e )/2  . X < 1/2 

' ''v  — — 


Equations  (D-30)  and  (D-31)  correspond  Co  Equations  (71)  and  (72)  of  Section 
2,6  These  equations  are  valid  only  in  the  left  half  of  the  vent  region 
((1  -P^)/2  < X < 1/2),  while  the  corresponding  values  In  the  right  half  of 

Region  3 are  obtained  by  symmetry. 
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